Small RPN exponential routine

[Gerson W. Barbosa]

This would make more sense on the HP-16C, which lacks exponentiation, but the equivalent program is 23 steps long.
I wrote a 42S version as an optimization exercise and managed to save five steps and two registers (but only two bytes), when compared to the first version. Perhaps more optimization is still possible, but I'll stop here.
I'll publish the HP-16C version, in case someone is interested (on the 16C this should be a first step towards a usable exponential function implementation).

Code:


00 { 29-Byte Prgm }
01>LBL "EX"
02 RCL ST Y
03 1
04>LBL 00
05 RCL× ST Z
06 STO+ ST Y
07 X<>Y
08 RCL× ST T
09 X<>Y
10 DSE ST Z
11 GTO 00
12 ÷
13 RCL÷ ST Z
14 .END.

Examples:

1 ENTER 14 XEQ EX --> 2.71828182846

3.14159265359 ENTER 22 --> 23.1406926328

1000 ENTER 1400 --> 1.97007111402e434 ( Free 42 only! )

[Thomas Klemm]

(01-02-2015 10:25 PM)Gerson W. Barbosa Wrote:  Perhaps more optimization is still possible, but I'll stop here.

Is there a reason for not using the obvious?

Code:

00 { 22-Byte Prgm }
01>LBL "EX"
02 1
03>LBL 00
04 RCL× ST Z
05 RCL÷ ST Y
06 1
07 +
08 DSE ST Y
09 GTO 00
10 END

Cheers
Thomas

[Gerson W. Barbosa]

(01-03-2015 08:42 AM)Thomas Klemm Wrote:  Is there a reason for not using the obvious?

Yes, dumbness...

D'oh!

Cheers,

Gerson.

[Gerson W. Barbosa]

This is my first attempt at natural logarithm using

ln(1 + x) = x - x^2/2 + x^3/3 - x^4/4 + x^5/5 -+ ..., x < 1

There surely is a more obvious method which I cannot see, so I won't bother optimizing it. I am primarily interested in implementing these functions on TurboBCD, but rather than trying to reinvent the wheel I'll look for a ready-made solution.

Code:


00 { 43-Byte Prgm }
01>LBL "LX"
02 NOT
03 +/-
04 LASTX
05 1
06 RCL- ST T
07 RCL÷ ST Z
08 LASTX
09 1
10>LBL 00
11 RCL÷ ST T
12 RCL+ ST Z
13 RCL× ST Y
14 STO ST Z
15 Rv
16 1
17 DSE ST T
18 GTO 00
19 RCL ST Z
20 +/-
21 .END.

0.52 ENTER 35 --> -6.53926467407e-1

There are two problems, though:

1) The series does not converge for arguments greater than 1;

2) More and more terms are required as arguments approach to zero. For instance, ln(0.052) will require 525 terms for 12-digit accuracy.

In these cases, conversions to and from decimal logarithms will help:


ln(0.52) = -6.53926467407e-1

2.30258509299 / --> -2.83996656366e-1

1 - 2.30258509299 * --> -2.95651156041 = ln(0.052)


ln(0.52) = -6.53926467407e-1

2.30258509299 / --> -2.83996656366e-1

1 + 2.30258509299 * --> 1.64865862558 = ln(5.2)


ln(0.52) = -6.53926467407e-1

2.30258509299 / --> -2.83996656366e-1

2 + 2.30258509299 * --> 3.95124371856 = ln(52)


ln(0.52) = -6.53926467407e-1

2.30258509299 / --> -2.83996656366e-1

3 + 2.30258509299 * --> 6.25382881155 = ln(520)

[Thomas Klemm]

(01-03-2015 05:31 PM)Gerson W. Barbosa Wrote:  There surely is a more obvious method which I cannot see, so I won't bother optimizing it.

Using the same method but a little shorter:

Code:

00 { 24-Byte Prgm }
01>LBL "LX"
02 1
03 STO- ST Z
04 CLX
05>LBL 00
06 RCL Y
07 1/X
08 X<>Y
09 -
10 RCL× ST Z
11 DSE ST Y
12 GTO 00
13 END

Cheers
Thomas

[Gerson W. Barbosa]

(01-03-2015 08:33 PM)Thomas Klemm Wrote:  

(01-03-2015 05:31 PM)Gerson W. Barbosa Wrote:  There surely is a more obvious method which I cannot see, so I won't bother optimizing it.

Using the same method but a little shorter:

Code:

00 { 24-Byte Prgm }
01>LBL "LX"
02 1
03 STO- ST Z
04 CLX
05>LBL 00
06 RCL Y
07 1/X
08 X<>Y
09 -
10 RCL× ST Z
11 DSE ST Y
12 GTO 00
13 END

Cheers
Thomas

A "little shorter"? I appreciate your sense of humor :-)

Definitely I should have made a second attempt in the first place. I've rewritten it looking at the expression with five terms of my example and obtained the same number of steps, but four bytes longer. Despite having one extra step inside the loop, yours appears to be a little faster (only one measurement before I got "Memory Clear" due to flat batteries).

Cheers,

Gerson.


ln(0.52) = -(0.48*(1 + 0.48*(1/2 + 0.48*(1/3 + 0.48*(1/4 + 0.48/5)))))

Code:


00 { 28-Byte Prgm }
01>LBL "LX"
02 1
03 RCL- ST Z
04 0
05>LBL 00
06 1
07 RCL÷ ST T
08 +
09 RCL× ST Y
10 DSE ST Z
11 GTO 00
12 +/-
13 .END.

P.S.: One step less, but still one byte longer, using one of your ideas:

Code:


00 { 25-Byte Prgm }
01>LBL "LX"
02 1
03 RCL- ST Z
04 0
05>LBL 00
06 RCL ST Z
07 1/X
08 -
09 RCL× ST Y
10 DSE ST Z
11 GTO 00
12 END

[James Dunn]

Gerson, Is your exponential routine similar to Toth's? (listed below)

Code:

001 - 43,22, d  LBL D
002 - 44 8      STO 8
003 - 1         1
004 - 44 32     STO I
005 - 44 9      STO 9
006 - 20        ×
007 - 43,22, 4  LBL 4
008 - 36        ENTER
009 - 36        ENTER
010 - 45 9      RCL 9
011 - 40        +
012 - 45 9      RCL 9
013 - 43 49     x=y
014 - 43 21     RTN
015 - 33        Rv
016 - 44 9      STO 9
017 - 43 24     ISZ
018 - 33        Rv
019 - 45 8      RCL 8
020 - 20        ×
021 - 45 32     RCL I
022 - 10        ÷
023 - 22 4      GTO 4

It works, but I wonder how. Could someone explain? Thanks.

[Namir]

How about HP-41C versions .... pleeeeeeeease!!

:-)

Namir

[Gerson W. Barbosa]

(01-20-2015 07:57 PM)James Dunn Wrote:  Gerson, Is your exponential routine similar to Toth's? (listed below)

Code:

001 - 43,22, d  LBL D
002 - 44 8      STO 8
003 - 1         1
004 - 44 32     STO I
005 - 44 9      STO 9
006 - 20        ×
007 - 43,22, 4  LBL 4
008 - 36        ENTER
009 - 36        ENTER
010 - 45 9      RCL 9
011 - 40        +
012 - 45 9      RCL 9
013 - 43 49     x=y
014 - 43 21     RTN
015 - 33        Rv
016 - 44 9      STO 9
017 - 43 24     ISZ
018 - 33        Rv
019 - 45 8      RCL 8
020 - 20        ×
021 - 45 32     RCL I
022 - 10        ÷
023 - 22 4      GTO 4

My HP-16C version:

Code:


01 g LBL A
02   STO I
03   1
04   STO 0
05 g CLx
06   x<>y
07   ENTER
08   ENTER
09   ENTER
10 g LBL 0
11   RCL 0
12   RCL I
13   *
14   STO 0
15   + 
16   *
17 g DSZ
20   GTO 0
19   RCL 0
20   /
21   x<>y
22   /
23 g RTN

Both have approximately the same size and both use Horner's method. My first HP-16C version of Thomas Klemm's more obvious method is only one step shorter, but I surely haven't converted it properly.

Regards,

Gerson.

P.S.: Contrary to what I stated, Victor Toth's program uses a different method. My bad, sorry!

[Namir]

(01-20-2015 11:48 PM)Gerson W. Barbosa Wrote:  My HP-16C version:

Code:


01 g LBL A
02   STO I
03   1
04   STO 0
05 g CLx
06   x<>y
07   ENTER
08   ENTER
09   ENTER
10 g LBL 0
11   RCL 0
12   RCL I
13   *
14   STO 0
15   + 
16   *
17 g DSZ
20   GTO 0
19   RCL 0
20   /
21   x<>y
22   /
23 g RTN

Both have approximately the same size and both use Horner's method. My first HP-16C version of Thomas Klemm's more obvious method is only one step shorter, but I surely haven't converted it properly.

Regards,

Gerson.

I think LBL A needs to be followed by an X<>Y command. I am assuming the value for the exponent is in register Y and the number of iterations is in register X. I tried your code on an HP-41C and it worked after inserting X<>Y following LBL A.

Namir

[Gerson W. Barbosa]

(01-21-2015 12:34 AM)Namir Wrote:  

(01-20-2015 11:48 PM)Gerson W. Barbosa Wrote:  My HP-16C version:

Code:


01 g LBL A
02   STO I
03   1
04   STO 0
05 g CLx
06   x<>y
07   ENTER
08   ENTER
09   ENTER
10 g LBL 0
11   RCL 0
12   RCL I
13   *
14   STO 0
15   + 
16   *
17 g DSZ
20   GTO 0
19   RCL 0
20   /
21   x<>y
22   /
23 g RTN

Both have approximately the same size and both use Horner's method. My first HP-16C version of Thomas Klemm's more obvious method is only one step shorter, but I surely haven't converted it properly.

Regards,

Gerson.

I think LBL A needs to be followed by an X<>Y command. I am assuming the value for the exponent is in register Y and the number of iterations is in register X. I tried your code on an HP-41C and it worked after inserting X<>Y following LBL A.

Namir

Namir,

On the 16C, 1 ENTER 12 GSB A returns 2.718281829. No x<>y need. I'll try it on the 41C.

Gerson.

P.S.: Instead of X<>Y after LBL A, try

STO 01 ; assuming that's your index register;
SIGN
STO 00
X<>Y
ENTER
ENTER
ENTER
...

Of course converting Thomas Klemm's program to the HP-41C is better.

[Thomas Klemm]

(01-20-2015 07:57 PM)James Dunn Wrote:  Gerson, Is your exponential routine similar to Toth's?
It works, but I wonder how. Could someone explain? Thanks.

Toth's program calculates the partial sums:

\(1\)
\(1+x\)
\(1+x+\frac{x^2}{2}\)
\(1+x+\frac{x^2}{2}+\frac{x^3}{6}\)
\(1+x+\frac{x^2}{2}+\frac{x^3}{6}+\frac{x^4}{24}\)
(...)

It returns if the value doesn't change anymore. This value is kept in register 9 while the original x is kept in register 8. The terms to be added are calculated recursively from the previous value: \(a_n=\frac{x}{n}a_{n-1}\) (lines 019-021). The index \(n\) is kept in register I.

Horner's method let us rewrite the expression:

\(1+x(1+\frac{x}{2}(1+\frac{x}{3}(1+\frac{x}{4}(...)))\)

My program calculates that inside out using classic RPN style. Thus the loop is inverted: the index is going down to 0. From a numerical point of view that's favorable as it reduces cancelation: the values that are added are of a similar size.

Best is to step through the program using special values like 2 that are easy to follow.

HTH
Thomas

[Thomas Klemm]

(01-20-2015 11:48 PM)Gerson W. Barbosa Wrote:  My first HP-16C version of Thomas Klemm's more obvious method is only one step shorter, but I surely haven't converted it properly.

This is how I would write my program for the HP-16C:

Code:

001 - 43,22, A  LBL A
002 - 44 32     STO I
003 - 33        X<>Y
004 - 44 0      STO 0
005 - 1         1
006 - 43,22, 0  LBL 0
007 - 45 0      RCL 0
008 - 20        ×
009 - 45 32     RCL I
010 - 10        ÷
011 - 1         1
012 - 40        +
013 - 43 23     DSZ
014 - 22 0      GTO 0
015 - 43 21     RTN

Cheers
Thomas

[Gerson W. Barbosa]

(01-21-2015 03:18 AM)Thomas Klemm Wrote:  

(01-20-2015 11:48 PM)Gerson W. Barbosa Wrote:  My first HP-16C version of Thomas Klemm's more obvious method is only one step shorter, but I surely haven't converted it properly.

This is how I would write my program for the HP-16C:

Code:

001 - 43,22, A  LBL A
002 - 44 32     STO I
003 - 33        X<>Y
004 - 44 0      STO 0
005 - 1         1
006 - 43,22, 0  LBL 0
007 - 45 0      RCL 0
008 - 20        ×
009 - 45 32     RCL I
010 - 10        ÷
011 - 1         1
012 - 40        +
013 - 43 23     DSZ
014 - 22 0      GTO 0
015 - 43 21     RTN

Cheers
Thomas

I'm glad my 22-step version has gone straight to the trash bin :-)

Cheers,

Gerson.

[Gerson W. Barbosa]

(01-20-2015 09:24 PM)Namir Wrote:  How about HP-41C versions .... pleeeeeeeease!!

:-)

Namir

Conversions from Thomas Klemm's HP-42S programs above are straightforward:

Code:


01 LBL 'LN                  01 LBL 'EX
02 1                        02 1            
03 ST- Z                    03 LBL 00
04 CLX                      04 RCL Z
05 LBL 00                   05 *
06 RCL Y                    06 RCL Y
07 1/X                      07 /
08 X<>Y                     08 1
09 -                        09 +
10 RCL Z                    10 DSE Y
11 *                        11 GTO 00
12 DSE Y                    12 END
13 GTO 00
14 END

Gerson.

[Namir]

Thank you for the HP-41C Code. Now for a slightly more difficult request. How about a version fo rthe HP-67?

Namir

[Gerson W. Barbosa]

(01-22-2015 08:41 AM)Namir Wrote:  Thank you for the HP-41C Code. Now for a slightly more difficult request. How about a version fo rthe HP-67?

Namir

No HP-67 here, but I have a printed manual. These should be equivalent to the above HP-41C programs and hopefully work on the HP-67:

Code:


001 h LBL A                    016 h LBL B
002 h ST I                     017 h ST I
003 h x<>y                     018 h x<>y
004   STO 1                    019   1
005   1                        020   -
006 h LBL 1                    021   STO 1
007   RCL 1                    022   0
008   *                        023 h LBL 2
009 h RC I                     024 h RC I
010   /                        025   1/x
011   1                        026 h x<>y
012   +                        027   -
013 f DSZ                      028   RCL 1
014   GTO 1                    029   *
015 h RTN                      030 f DSZ
                               031   GTO 2
                               032 h RTN

Gerson.

[Namir]

(01-22-2015 06:59 PM)Gerson W. Barbosa Wrote:  

(01-22-2015 08:41 AM)Namir Wrote:  Thank you for the HP-41C Code. Now for a slightly more difficult request. How about a version fo rthe HP-67?

Namir

No HP-67 here, but I have a printed manual. These should be equivalent to the above HP-41C programs and hopefully work on the HP-67:

Code:


001 h LBL A                    016 h LBL B
002 h ST I                     017 h ST I
003 h x<>y                     018 h x<>y
004   STO 1                    019   1
005   1                        020   -
006 h LBL 1                    021   STO 1
007   RCL 1                    022   0
008   *                        023 h LBL 2
009 h RC I                     024 h RC I
010   /                        025   1/x
011   1                        026 h x<>y
012   +                        027   -
013 f DSZ                      028   RCL 1
014   GTO 1                    029   *
015 h RTN                      030 f DSZ
                               031   GTO 2
                               032 h RTN

Gerson.

Gerson,

Thank you so much for the HP-67 code.

Have you considered attending an HHC conference in the US? It will be a special treat to meet you in person. I am sure you will have a time of your life being among peers who respect you and your contributions here.

Namir

[Gerson W. Barbosa]

(01-23-2015 05:11 PM)Namir Wrote:  Gerson,

Thank you so much for the HP-67 code.

Have you considered attending an HHC conference in the US? It will be a special treat to meet you in person. I am sure you will have a time of your life being among peers who respect you and your contributions here.

Namir

Namir,

Thank you very much for your kind words, despite my ever fading contributions.
Yes, I've already considered attendind one HHC conference. It would be a pleasure meeting you all in person. Perhaps one of these years :-)

Thanks again for the invitation!

Gerson.

[Eddie W. Shore]

Is this the origin for the LPN1 function? (I need to be here more often)