Natural logarithm of 2 [HP-15C/HP-42S/Free42 & others]

[Gerson W. Barbosa]

My goal is finding a fast and short way to compute the natural logarithm of 2 on the HP-15C and other RPN calculators like the HP-42S, without using any built-in function (LN, LOG, ATANH...). By fast I mean under 10 seconds on the HP-15C for maximum possible precision. By short I mean no more than 30 or 40 RPN steps.

I have provided two RPN programs and a Decimal BASIC program. The latter avoids any multiplication operations inside the loop, except for the two necessary divisions. I would like to have tried it on the HP-11C which is slightly faster, but currently I have none.

Later I will post the (partially original) formula I have used, in case it is not so obvious from the listings.

You are invited to contribute code for other RPN calculators (HP-41C, HP-32S, HP-12C...), based on this algorithm or improved versions thereof, or based on more efficient ones (plenty of formulas in the above hyperlink). Pseudo code and short BASIC programs are also welcome.

Thanks in advance for your interest and contribuitions.

Gerson.

======================================

HP-15C

======================================


001-    42 21 11  f LBL A
002-       44 25    STO I
003-           4    4
004-          20    ×
005-           1    1
006-          40    +
007-       44  1    STO 1
008-           0    0
009-       44  0    STO 0
010-    42 21  0  f LBL 0
011-    45 40  1    RCL + 1
012-       45 25    RCL I
013-       43 11  g x²
014-          34    x↔y
015-          10    ÷
016-       45 25    RCL I
017-    45 40 25    RCL + I
018-       43 11  g x²
019-       43 36  g LSTx
020-          30    -
021-          15    1/x
022-    44 40  0    STO + 0
023-          34    x↔y
024-    42  5 25  f DSE I
025-       22  0    GTO 0
026-    45 40  1    RCL + 1
027-          15    1/x
028-    45 40  0    RCL + 0
029-       43 32  g RTN



4 f A -> 0.6931471808 ( ~  8.5 seconds )

5 f A -> 0.6931471805 ( ~ 10.0 seconds )


======================================

HP-42S/Free42

======================================


00 { 48-Byte Prgm }
01▸LBL "Ln2"
02 RCL ST X
03 4
04 ×
05 1
06 +
07 0
08 STO 00
09▸LBL 00
10 RCL+ ST Y
11 RCL ST Z
12 X↑2
13 X<>Y
14 ÷
15 RCL ST Z
16 STO+ ST X
17 X↑2
18 RCL- ST L
19 1/X
20 STO+ 00
21 R↓
22 DSE ST Z
23 GTO 00
24 +
25 1/X
26 RCL+ 00
27 END

 6 XEQ Ln2 -> 0.69314718056 ( ~ 2 seconds, HP-42S ) 

16 XEQ Ln2 -> 0.6931471805599453094172321214581766 ( Free42 )

======================================

Decimal BASIC

======================================

OPTION ARITHMETIC DECIMAL_HIGH
! Log(2)
INPUT PROMPT "Number of decimal places: ":nd
LET t = TIME 
LET k = CEIL(ND/LOG(8))                                                         
LET n = k*k
LET m = 2*k - 1                                                   
LET c = 8                                                                                                                                                                                             
LET d0 = c*(k - 2) + 10 
LET d1 = 2*k*(2*k - 1)
LET d2 = 4*k + 1                                         
LET s1 = 0  
LET s2 = 0                                                                                                 
FOR i = 1 TO k                                                                                       
   LET s1 = s1 + 1/d1                                                  
   LET d1 = d1 - d0                     
   LET d0 = d0 - c
   LET s2 = n/(s2 + d2)                                               
   LET n = n - m
   LET m = m - 2
NEXT i
LET s2 = 1/(s2 + d2)
LET r = TIME - t 
LET r$ = STR$(s1 + s2)                           
PRINT
PRINT " ";
PRINT r$(0:1);
FOR i = 2 TO nd + 1
   PRINT r$(i:i);
   IF MOD((i - 1),10) = 0 THEN PRINT " ";
   IF MOD((i - 1),50) = 0 THEN 
      PRINT
      PRINT "  ";
   END IF
NEXT i
IF MOD (i - 2,50) <> 0  OR nd = 0 THEN PRINT
PRINT 
PRINT "Runtime: ";
PRINT  USING "0.##": r;
PRINT " seconds"
END


Number of decimal places: 1000

 .6931471805 5994530941 7232121458 1765680755 0013436025 
  5254120680 0094933936 2196969471 5605863326 9964186875 
  4200148102 0570685733 6855202357 5813055703 2670751635 
  0759619307 2757082837 1435190307 0386238916 7347112335 
  0115364497 9552391204 7517268157 4932065155 5247341395 
  2588295045 3007095326 3666426541 0423915781 4952043740 
  4303855008 0194417064 1671518644 7128399681 7178454695 
  7026271631 0645461502 5720740248 1637773389 6385506952 
  6066834113 7273873722 9289564935 4702576265 2098859693 
  2019650585 5476470330 6793654432 5476327449 5125040606 
  9438147104 6899465062 2016772042 4524529612 6879465461 
  9316517468 1392672504 1038025462 5965686914 4192871608 
  2938031727 1436778265 4877566485 0856740776 4845146443 
  9940461422 6031930967 3540257444 6070308096 0850474866 
  3852313818 1676751438 6674766478 9088143714 1985494231 
  5199735488 0375165861 2753529166 1000710535 5824987941 
  4729509293 1138971559 9820565439 2871700072 1808576102 
  5236889213 2449713893 2037843935 3088774825 9701715591 
  0708823683 6275898425 8918535302 4363421436 7061189236 
  7891923723 1467232172 0534016492 5687274778 2344535347 
  
Runtime: 0.32 seconds

[JMBaillard]

Hi !
Here is an HP-41 version:

01 LBL "LN2"
02 1
03 STO 00
04 2
05 STO 01
06 1/X
07 LBL 01
08 RCL 01
09 RCL 00
10 ST/ Y
11 1
12 +
13 STO 00
14 *
15 ST+ X
16 STO 01
17 1/X
18 +
19 X#Y?
20 GTO 01
21 END

(32 bytes)

It takes about 13 seconds with an HP-41 which gives 0.6931471808
So, an error of about 2 E-10 ( and perhaps a little too slow )
But with free42, the result 6,931471805599453094172321214581773e-1
is returned in much less than 1 second.

Best regards.

[John Keith]

(10-12-2019 12:32 AM)Gerson W. Barbosa Wrote:  Number of decimal places: 1000

 .6931471805 5994530941 7232121458 1765680755 0013436025 
  5254120680 0094933936 2196969471 5605863326 9964186875 
  4200148102 0570685733 6855202357 5813055703 2670751635 
  0759619307 2757082837 1435190307 0386238916 7347112335 
  0115364497 9552391204 7517268157 4932065155 5247341395 
  2588295045 3007095326 3666426541 0423915781 4952043740 
  4303855008 0194417064 1671518644 7128399681 7178454695 
  7026271631 0645461502 5720740248 1637773389 6385506952 
  6066834113 7273873722 9289564935 4702576265 2098859693 
  2019650585 5476470330 6793654432 5476327449 5125040606 
  9438147104 6899465062 2016772042 4524529612 6879465461 
  9316517468 1392672504 1038025462 5965686914 4192871608 
  2938031727 1436778265 4877566485 0856740776 4845146443 
  9940461422 6031930967 3540257444 6070308096 0850474866 
  3852313818 1676751438 6674766478 9088143714 1985494231 
  5199735488 0375165861 2753529166 1000710535 5824987941 
  4729509293 1138971559 9820565439 2871700072 1808576102 
  5236889213 2449713893 2037843935 3088774825 9701715591 
  0708823683 6275898425 8918535302 4363421436 7061189236 
  7891923723 1467232172 0534016492 5687274778 2344535347 
  
Runtime: 0.32 seconds

For comparison, the HP 50g with LongFloat returns the result above (except for the last digit being 8) in 637 seconds.

[Paul Dale]

The WP 34S algorithm for logarithms might be appropriate here:

Code:

/* Natural logarithm.
 *
 * Take advantage of the fact that we store our numbers in the form: m * 10^e
 * so log(m * 10^e) = log(m) + e * log(10)
 * do this so that m is always in the range 0.1 <= m < 2.  However if the number
 * is already in the range 0.5 .. 1.5, this step is skipped.
 *
 * Then use the fact that ln(x^2) = 2 * ln(x) to range reduce the mantissa
 * into 1/sqrt(2) <= m < 2.
 *
 * Finally, apply the series expansion:
 *   ln(x) = 2(a+a^3/3+a^5/5+...) where a=(x-1)/(x+1)
 * which converges quickly for an argument near unity.
 */

Either use an argument of 2 directly or take advantage of \( \frac{\sqrt{2} - 1}{\sqrt{2} + 1} = 3 - 2 \sqrt{2} \).


Pauli

[ttw]

Perhaps the methods used here may be of interest.

http://elib.mi.sanu.ac.rs/files/journals...tm1212.pdf

[Gerson W. Barbosa]

(10-12-2019 11:04 AM)JMBaillard Wrote:  Hi !
Here is an HP-41 version:

01 LBL "LN2"
02 1
03 STO 00
04 2
05 STO 01
06 1/X
07 LBL 01
08 RCL 01
09 RCL 00
10 ST/ Y
11 1
12 +
13 STO 00
14 *
15 ST+ X
16 STO 01
17 1/X
18 +
19 X#Y?
20 GTO 01
21 END

(32 bytes)

It takes about 13 seconds with an HP-41 which gives 0.6931471808
So, an error of about 2 E-10 ( and perhaps a little too slow )
But with free42, the result 6,931471805599453094172321214581773e-1
is returned in much less than 1 second.

Hello, Jean-Marc,

Very concise program!

Here is another one, based on formula #20 at MathWorld (not what I’ve been using).

01 LBL "LN2"
02 0
03 3
04 STO 00
05 9
06 STO 01
07 LBL 00
08 *
09 1/X
10 +
11 2
12 ST+ 00
13 RDN
14 9
15 ST* 01
16 RDN
17 RCL 00
18 RCL 01
19 DSE T
20 GTO 00
21 2
22 ST* T
23 R^
24 +
25 3
26 /
27 END

9 XEQ LN2 -> 0.6931471807 (6.5 seconds on my HP-41C)

I am sure this can be optimized for size, but the lack of recall arithmetic is always a problem to me.

Best regards,

Gerson.

[JMBaillard]

Hello Gerson,

thank you for the link, I didn't know this forrmula !
Here is another HP41-version:

01 LBL "LN2"
02 CLX
03 STO 00
04 SIGN
05 STO 01
06 LBL 01
07 RCL 01
08 9
09 *
10 STO 01
11 RCL 00
12 1
13 +
14 STO 00
15 ST+ X
16 LASTX
17 +
18 *
19 1/X
20 +
21 X#Y?
22 GTO 01
23 ST+ X
24 3
25 /
26 END

(37 bytes)

The HP41 returns 0.6931471807 in 5.1 seconds
and free42 gives 0,693147180559945309417232121458174

Best regards,
JM.

[Gerson W. Barbosa]

Here is an HP-41 equivalent of my previous HP-15C and HP-42S programs:


01 LBL "LN2"
02 RCL X
03 4
04 *
05 1
06 +
07 0
08 STO 00
09 LBL 00
10 RCL Y
11 +
12 RCL Z
13 X^2
14 X<>Y
15 /
16 R^
17 ST+ T
18 X<> T
19 X^2
20 ST- L
21 X<> L
22 1/X
23 ST- 00
24 RDN
25 DSE Z
26 GTO 00
27 +
28 1/X
29 RCL 00
30 +
31 END


4 XEQ LN2 -> 0.6931471808 ( ~ 2.9 seconds )

5 XEQ LN2 -> 0.6931471805 ( ~ 3.5 seconds )


P.S.: For the sake of completeteness, here is the formula I've been using:

\[\ln (2)\approx \sum_{n=1}^{k}\frac{1}{2n\left ( n+1 \right )}+\frac{1}{4k+1+\frac{1}{4k+1+\frac{4}{4k+1+\frac{9}{4k+1+\frac{\ddots }{ 4k+1+\frac{k^{2}}{4k+1}}}}}}\]

The known series together with this continued fraction correction terms yield about 2.09 digits per k. This comes from sheer observation, thus being provided with no proof. However, it appears to hold.

[toml_12953]

Here's a Python program modeled after the Decimal BASIC version but it doesn't have extended precision. I'm new to Python and don't know how to get the most out of it.

Code:

def log_2():
    import math,time
    nd=int(input("Number of decimal places: "))
    t = time.time()
    k = math.ceil(nd/math.log(8))
    n = k*k
    m = 2*k - 1
    c = 8
    d0 = c*(k - 2) + 10
    d1 = 2*k*(2*k - 1)
    d2 = 4*k + 1
    s1 = 0
    s2 = 0
    for i in range(1,k+1):
        s1 = s1 + 1/d1
        d1 = d1 - d0
        d0 = d0 - c
        s2 = n/(s2 + d2)
        n = n - m
        m = m - 2
    s2 = 1/(s2 + d2)
    r = time.time() - t
    rs = str(s1 + s2)
    print()
    print(" ",end="")
    print(rs[0:2],end="")
    for i in range(2,nd+2):
        print(rs[i:i+1],end="")
        if (i - 1) % 10 == 0:
            print(" ",end="")
        if (i - 1) % 50 == 0:
            print()
            print(" ",end="")
    if (i - 2) % 50 != 0 or nd == 0:
        print()
    print()
    print("Runtime: ",end="")
    print('%4.2f' %r ,end="")
    print(" seconds")

[Gerson W. Barbosa]

(10-12-2019 10:22 PM)toml_12953 Wrote:  Here's a Python program modeled after the Decimal BASIC version but it doesn't have extended precision. I'm new to Python and don't know how to get the most out of it.

Code:

def log_2():
    import math,time
    nd=input("Number of decimal places: ")
    nd=int(nd)
    t = time.time()
    k = math.ceil(nd/math.log(8))
    n = k*k
    m = 2*k - 1
    c = 8
    d0 = c*(k - 2) + 10
    d1 = 2*k*(2*k - 1)
    d2 = 4*k + 1
    s1 = 0
    s2 = 0
    for i in range(1,k+1):
        s1 = s1 + 1/d1
        d1 = d1 - d0
        d0 = d0 - c
        s2 = n/(s2 + d2)
        n = n - m
        m = m - 2
    s2 = 1/(s2 + d2)
    r = time.time() - t
    rs = str(s1 + s2)
    print()
    print(" ",end="")
    print(rs[0:2],end="")
    for i in range(2,nd+2):
        print(rs[i:i+1],end="")
        if (i - 1) % 10 == 0:
            print(" ",end="")
        if (i - 1) % 50 == 0:
            print()
            print(" ",end="")
    if (i - 2) % 50 != 0 or nd == 0:
        print()
    print()
    print("Runtime: ",end="")
    print('%4.2f' %r ,end="")
    print(" seconds")

Decimal BASIC is easy to use, but it’s limited to only 1000 digits.
I’ve noticed this algorithm produces about 2.09 digits per iteration, hence the use of ln(8) in the estimation of the required number of iterations for a given number of digits.

[Albert Chan]

Here is HP-71B code to calculate log(2)

Line 10: the code actually do log(8), with guess = 2
Line 20: Y = SINH(X-2) ≈ (X-2) + (X-2)^3/3! + (X-2)^5/5! + (X-2)^7/7!
Line 30: Y = EXP(2) * EXP(X-2) = EXP(X)
Line 40: cubically convergent iterations for LOG(N)

Code:

10 N=8 @ X=2 @ E2=EXP(2)
20 Y=X-2 @ Z=Y*Y @ Z=Z/6*(1+Z/20*(1+Z/42)) @ Y=Y+Y*Z
30 Y=E2*(Y+SQRT(1+Y*Y))
40 Z=2*(Y-N)/(Y+N) @ X=X-Z
45 PRINT X/3
50 IF ABS(Z)>.0001 THEN 20

>RUN
.69313326289
.69314718056

Ref: http://rnc7.loria.fr/brent_invited.pdf, page 20
Ref: https://pdfs.semanticscholar.org/2db2/10...5e12fe.pdf

[Gerson W. Barbosa]

(10-12-2019 12:22 PM)John Keith Wrote:  

(10-12-2019 12:32 AM)Gerson W. Barbosa Wrote:  Number of decimal places: 1000

 .6931471805 5994530941 7232121458 1765680755 0013436025 
  5254120680 0094933936 2196969471 5605863326 9964186875 
  4200148102 0570685733 6855202357 5813055703 2670751635 
  0759619307 2757082837 1435190307 0386238916 7347112335 
  0115364497 9552391204 7517268157 4932065155 5247341395 
  2588295045 3007095326 3666426541 0423915781 4952043740 
  4303855008 0194417064 1671518644 7128399681 7178454695 
  7026271631 0645461502 5720740248 1637773389 6385506952 
  6066834113 7273873722 9289564935 4702576265 2098859693 
  2019650585 5476470330 6793654432 5476327449 5125040606 
  9438147104 6899465062 2016772042 4524529612 6879465461 
  9316517468 1392672504 1038025462 5965686914 4192871608 
  2938031727 1436778265 4877566485 0856740776 4845146443 
  9940461422 6031930967 3540257444 6070308096 0850474866 
  3852313818 1676751438 6674766478 9088143714 1985494231 
  5199735488 0375165861 2753529166 1000710535 5824987941 
  4729509293 1138971559 9820565439 2871700072 1808576102 
  5236889213 2449713893 2037843935 3088774825 9701715591 
  0708823683 6275898425 8918535302 4363421436 7061189236 
  7891923723 1467232172 0534016492 5687274778 2344535347 
  
Runtime: 0.32 seconds

For comparison, the HP 50g with LongFloat returns the result above (except for the last digit being 8) in 637 seconds.

Considering 0.32 seconds was achieved on a 2.6 GHz, that’s pretty good on the HP-50g @ 125 Mhz or so and one emulation layer.

[Paul Dale]

(10-12-2019 04:18 PM)Gerson W. Barbosa Wrote:  Here is another one, based on formula #20 at MathWorld (not what I’ve been using).

This looks to be the same as the method the 34S uses but calculating the series at x=2 directly (i.e. without the prescaling).


Pauli

[Gerson W. Barbosa]

HP 50g:

Code:


%%HP: T(3)A(R)F(.);
\<< RCLF SWAP RAD -105 CF -3 CF R\->I DUP 8 LN / CEIL 8 OVER 2 - OVER * 10 + PICK3 DUP + SQ LASTARG - 4 PICK DUP
  SQ OVER DUP + 1 - ROT 4 * 1 + 0 DUP \-> c d0 d1 n m d2 s1 s2
  \<< 1 SWAP
    START d1 INV 's1' STO+ d0 NEG 'd1' STO+ c NEG 'd0' STO+ n s2 d2 + / 's2' STO m NEG 'n' STO+ -2 'm' STO+
    NEXT s1 d2 s2 +
  \>> INV + EXPAND FXND ROT ALOG OVER - PICK3 * SWAP IQUOT + \->STR -51 FC? { "0." } { "0," } IFTE SWAP + SWAP STOF
\>>

200 -> 0.69314718055994530941723212145817656807550013436025
         52541206800094933936219696947156058633269964186875
         42001481020570685733685520235758130557032670751635
         07596193072757082837143519030703862389167347112335

( 350 seconds )

[Albert Chan]

Using Pi and AGM to get log(2): see http://rnc7.loria.fr/brent_invited.pdf, page 32

\(\large log(2^n) = \frac{\pi/2}{AGM(1, \frac{4}{2^n})} \left(1 + O(\frac{1}{2^{2n}}) \right)\)

Big enough n so that error term may be ignored. Assume O(1/4^n) = 1/4^n :

1/4^n < 1e-12
n ≥ 12 / log10(4) ≥ 19.93

Casio FX-115MS, with n = 20:

A=1
B=2^-18
C=A+B : B=√AB : A=C/2 : A-B

====     → 0.4980487824
====     → 0.2197274566
====     → 0.0525527223
====     → 0.0030466068
====     → 0.0000102395
====     → 0.0000000001

Convergence close to quadratic. With final A,B, we have

log(2) ≈ pi/(n(A+B)) = 0.693147180558, under-estimated 2e-12

Update: More precise calculations shows that n=20 over-estimated log(2) by 2e-12
Below suggested relative error term O(1/4^n) ≈ 1/4^(n-1)

python> from gmpy2 import *
python> f = lambda n: const_pi() / (2 * n * agm(1, 4/2**n))
python> for n in range(20, 25): print n, "%.15g" % f(n)
...
20 0.693147180562285
21 0.693147180560532
22 0.693147180560093
23 0.693147180559982
24 0.693147180559955

[Gerson W. Barbosa]

(11-08-2019 03:18 PM)Albert Chan Wrote:  Using Pi and AGM to get log(2): see http://rnc7.loria.fr/brent_invited.pdf, page 32

\(\large log(2^n) = \frac{\pi/2}{AGM(1, \frac{4}{2^n})} \left(1 + O(\frac{1}{2^{2n}}) \right)\)

Big enough n so that error term may be ignored. Assume O(1/4^n) = 1/4^n :

1/4^n < 1e-12
n ≥ 12 / log10(4) ≥ 19.93

Casio FX-115MS, with n = 20:

A=1
B=2^-18
C=A+B : B=√AB : A=C/2 : A-B

====     → 0.4980487824
====     → 0.2197274566
====     → 0.0525527223
====     → 0.0030466068
====     → 0.0000102395
====     → 0.0000000001

Convergence close to quadratic. With final A,B, we have

log(2) ≈ pi/(n(A+B)) = 0.693147180558, under-estimated 2e-12

Update: More precise calculations shows that n=20 over-estimated log(2) by 2e-12
Below suggested relative error term O(1/4^n) ≈ 1/4^(n-1)

python> from gmpy2 import *
python> f = lambda n: const_pi() / (2 * n * agm(1, 4/2**n))
python> for n in range(20, 25): print n, "%.15g" % f(n)
...
20 0.693147180562285
21 0.693147180560532
22 0.693147180560093
23 0.693147180559982
24 0.693147180559955

Interesting, even though it requires \(\pi\).

Because I’m lazy I’ll take the wp34s for this one:


001 LBL A
002 # π
003 RCL Y
004 # 002
005 -
006 # 1/2
007 STO+ X
008 DSE Y
009 BACK 002
010 STO* Z
011 STO+ X
012 INC Y
013 AGM
014 RCL* T
015 /
016 END

20 A -> 0.6931471805622850
2 LN - -> 2.33975581274e-12

40 A 2 LN - -> 2.210713576e-24

62 A 2 LN - -> 0


2^x would’ve saved a few steps, but I want to avoid the use of logarithms in the program.

[Gerson W. Barbosa]

(11-08-2019 03:18 PM)Albert Chan Wrote:  Update: More precise calculations shows that n=20 over-estimated log(2) by 2e-12

n = 1663 will be enough for 1000 digits:

--------------------------------------------------------

OPTION ARITHMETIC DECIMAL_HIGH
LET n = 1663
LET nd = 1000
PRINT "Log(2) = "
LET t = TIME 
LET p = 2^(n - 3)
LET d = PI*p
LET a = 1
LET b = p + p
FOR i = 1 TO 20
   LET c = a
   LET a = (a + b)/2
   LET b = SQR(b*c)
NEXT i
LET g = d/(a*n)
LET r = TIME - t
LET r$ = STR$(g)                           
PRINT
PRINT "0";
PRINT r$(0:1);
FOR i = 2 TO nd + 1
   PRINT r$(i:i);
   IF MOD((i - 1),10) = 0 THEN PRINT " ";
   IF MOD((i - 1),50) = 0 THEN 
      PRINT
      PRINT "  ";
   END IF
NEXT i
IF MOD (i - 2,50) <> 0  OR nd = 0 THEN PRINT
PRINT 
PRINT "Runtime: ";
PRINT  USING "0.##": r;
PRINT " seconds"
END

--------------------------------------------------------

Log(2) = 

0.6931471805 5994530941 7232121458 1765680755 0013436025 
  5254120680 0094933936 2196969471 5605863326 9964186875 
  4200148102 0570685733 6855202357 5813055703 2670751635 
  0759619307 2757082837 1435190307 0386238916 7347112335 
  0115364497 9552391204 7517268157 4932065155 5247341395 
  2588295045 3007095326 3666426541 0423915781 4952043740 
  4303855008 0194417064 1671518644 7128399681 7178454695 
  7026271631 0645461502 5720740248 1637773389 6385506952 
  6066834113 7273873722 9289564935 4702576265 2098859693 
  2019650585 5476470330 6793654432 5476327449 5125040606 
  9438147104 6899465062 2016772042 4524529612 6879465461 
  9316517468 1392672504 1038025462 5965686914 4192871608 
  2938031727 1436778265 4877566485 0856740776 4845146443 
  9940461422 6031930967 3540257444 6070308096 0850474866 
  3852313818 1676751438 6674766478 9088143714 1985494231 
  5199735488 0375165861 2753529166 1000710535 5824987941 
  4729509293 1138971559 9820565439 2871700072 1808576102 
  5236889213 2449713893 2037843935 3088774825 9701715591 
  0708823683 6275898425 8918535302 4363421436 7061189236 
  7891923723 1467232172 0534016492 5687274778 2344535347 
  
Runtime: 0.08 seconds

--------------------------------------------------------



That's one fourth of the time required by my original program based on series combined with a continued fraction. Only 20 iterations here compared to 481 there, but AGM requires the evaluation of one square root per iteration. I'll post later another solution based on a previously suggested algorithm which does it in 0.04 seconds. It would be better to compare those running times on a compiled language, though.

[Albert Chan]

(11-11-2019 06:14 PM)Gerson W. Barbosa Wrote:  Only 20 iterations here compared to 481 there, but AGM requires the evaluation of one square root per iteration.

We can lower loop count by 2, and extrapolate for the AGM.
(it probably make no difference to the timings though ...)

Loop 19 arithmetic mean: c = (a+b)/2
Loop 19 geometric mean: d = √(ab) = √(c^2 - ε^2), where ε = (a-b)/2

Taylor series for √(1+x) = 1 + x/2 - x^2/8 + x^3/16 - 5x^4/128 + ...

→ d = c * √(1 - (ε/c)^2) ≈ c - ε^2 / (2c)

Loop 20 aritmetic mean:

AGM ≈ ½(c + d) = c - ε^2 / (4c) = c - (a-b)^2 / (16*c)

Code:

FOR i = 1 TO 18
   LET c = a
   LET a = (a + b)/2
   LET b = SQR(b*c)
NEXT i
LET c = (a + b)/2
LET a = c - (a - b)^2 / (16*c)
LET g = d/(a*n)
...

[Gerson W. Barbosa]

(11-12-2019 08:30 PM)Albert Chan Wrote:  

(11-11-2019 06:14 PM)Gerson W. Barbosa Wrote:  Only 20 iterations here compared to 481 there, but AGM requires the evaluation of one square root per iteration.

We can lower loop count by 2, and extrapolate for the AGM.
(it probably make no difference to the timings though ...)

Loop 19 arithmetic mean: c = (a+b)/2
Loop 19 geometric mean: d = √(ab) = √(c^2 - ε^2), where ε = (a-b)/2

Taylor series for √(1+x) = 1 + x/2 - x^2/8 + x^3/16 - 5x^4/128 + ...

→ d = c * √(1 - (ε/c)^2) ≈ c - ε^2 / (2c)

Loop 20 aritmetic mean:

AGM ≈ ½(c + d) = c - ε^2 / (4c) = c - (a-b)^2 / (16*c)

Code:

FOR i = 1 TO 18
   LET c = a
   LET a = (a + b)/2
   LET b = SQR(b*c)
NEXT i
LET c = (a + b)/2
LET a = c - (a - b)^2 / (16*c)
LET g = d/(a*n)
...

Albert Chan, you rock!

OPTION ARITHMETIC DECIMAL_HIGH
LET n = 1663
LET nd = 1000
PRINT "Log(2) = "
LET t = TIME 
LET p = 2^(n - 3)
LET d = PI*p
LET a = 1
LET b = p + p
FOR i = 1 TO 18
   LET c = a
   LET a = (a + b)/2
   LET b = SQR(b*c)
NEXT i
LET c = (a + b)/2
LET a = c - (a - b)^2/(16*c)
LET g = d/(a*n)
LET r = TIME - t
LET r$ = STR$(g)                           
PRINT
PRINT "0";
PRINT r$(0:1);
FOR i = 2 TO nd + 1
   PRINT r$(i:i);
   IF MOD((i - 1),10) = 0 THEN PRINT " ";
   IF MOD((i - 1),50) = 0 THEN 
      PRINT
      PRINT "  ";
   END IF
NEXT i
IF MOD (i - 2,50) <> 0  OR nd = 0 THEN PRINT
PRINT 
PRINT "Runtime: ";
PRINT  USING "0.##": r;
PRINT " seconds"
END

----------------------------------

Log(2) = 

0.6931471805 5994530941 7232121458 1765680755 0013436025 
  5254120680 0094933936 2196969471 5605863326 9964186875 
  4200148102 0570685733 6855202357 5813055703 2670751635 
  0759619307 2757082837 1435190307 0386238916 7347112335 
  0115364497 9552391204 7517268157 4932065155 5247341395 
  2588295045 3007095326 3666426541 0423915781 4952043740 
  4303855008 0194417064 1671518644 7128399681 7178454695 
  7026271631 0645461502 5720740248 1637773389 6385506952 
  6066834113 7273873722 9289564935 4702576265 2098859693 
  2019650585 5476470330 6793654432 5476327449 5125040606 
  9438147104 6899465062 2016772042 4524529612 6879465461 
  9316517468 1392672504 1038025462 5965686914 4192871608 
  2938031727 1436778265 4877566485 0856740776 4845146443 
  9940461422 6031930967 3540257444 6070308096 0850474866 
  3852313818 1676751438 6674766478 9088143714 1985494231 
  5199735488 0375165861 2753529166 1000710535 5824987941 
  4729509293 1138971559 9820565439 2871700072 1808576102 
  5236889213 2449713893 2037843935 3088774825 9701715591 
  0708823683 6275898425 8918535302 4363421436 7061189236 
  7891923723 1467232172 0534016492 5687274778 2344535347 
  
Runtime: 0.04 seconds

My system:

Processor: Intel(R) Core (TM)2 Duo CPU E4700 @ 2.6GHz 2.6 GHz
64-bit Operational System

[Gerson W. Barbosa]

Former forum member Mike (Stgt) sent me two RPN programs which give perfect 10-digit and 12-digit results instantaneously on RPN calculators, inspired in Paul Dale's suggestion in message #4 of this thread:


01 LBL "LN2"
02 66
03 2
04 RCL Y
05 1/X
06 Y^X
07 RCL X
08 1
09 ST- Z
10 +
11 /
12 ENTER 
13 X^2
14 3
15 ST+ Y
16 /
17 *
18 *
19 ST+ X
20 END

That's log(2) = 132*atanh((2^(1/66) - 1)/(1 + 2^(1/66))), evaluated to two terms of the atanh(x) series:

132*((2^(1/66) - 1)/(1 + 2^(1/66)) + 1/3*((2^(1/66) - 1)/(1 + 2^(1/66)))^3) = 0.6931471806, on the HP-41


01▶LBL "LN2"
02 74
03 2
04 RCL ST Y
05 1/X
06 Y^X
07 RCL ST X
08 1
09 STO- ST Z
10 +
11 ÷
12 ENTER
13 X^2
14 3
15 STO+ ST Y
16 ÷
17 ×
18 ×
19 STO+ ST X
20 END

That's log(2) = 148*atanh((2^(1/74) - 1)/(1 + 2^(1/74))), evaluated to two terms of the atanh(x) series:

148*((2^(1/74) - 1)/(1 + 2^(1/74)) + 1/3*((2^(1/74) - 1)/(1 + 2^(1/74)))^3) = 0.693147180560, on the HP-42S

Inspired in Mike's contribution, I've written two HP-42S programs that give perfect results, either on the HP-42S or on Free42S:

---------------------------------------------------

00 { 52-Byte Prgm }
01▸LBL "Ln2"
02 3
03 2
04 SQRT
05 FP
06 1
07 RCL+ ST L
08 ÷
09 STO 00
10 X↑2
11 LASTX
12▸LBL 00
13 RCL× ST Y
14 RCL÷ ST Z
15 STO+ 00
16 X<> ST L
17 X<> ST Z
18 SIGN
19 STO+ ST X
20 RCL+ ST L
21 X<> ST Z
22 DSE ST T
23 GTO 00
24 4
25 RCL× 00
26 .END.

 6 XEQ Ln2 -> 0,69314718056 (under 2 seconds on the HP-42S)

21 XEQ Ln2 -> 0,6931471805599453094172321214581760 (Free42)
       
      2 LN -> 0,6931471805599453094172321214581766 (Free42)

---------------------------------------------------

00 { 55-Byte Prgm }
01▸LBL "Ln2"
02 3
03 2
04 SQRT
05 SQRT
06 1
07 -
08 2
09 RCL+ ST Y
10 ÷
11 STO 00
12 X↑2
13 LASTX
14▸LBL 00
15 RCL× ST Y
16 RCL÷ ST Z
17 STO+ 00
18 X<> ST L
19 X<> ST Z
20 SIGN
21 STO+ ST X
22 RCL+ ST L
23 X<> ST Z
24 DSE ST T
25 GTO 00
26 8               ; double this constant after each additional SQRT
27 RCL× 00
28 .END.


 4 XEQ Ln2 -> 0,69314718055 (HP-42S)

14 XEQ Ln2 -> 0,6931471805599453094172321214581766 (Free42)
       
      2 LN -> 0,6931471805599453094172321214581766 (Free42)

--------------------------------------------------- 

This approach has allowed me to get 999 correct digits in 0.03s on my computer, with Decimal BASIC:

 
OPTION ARITHMETIC DECIMAL_HIGH
LET nd = 1000
PRINT "Log(2) ="
LET t = TIME 
LET k = 298                                                   
LET r = TIME - t 
LET a = SQR(SQR(SQR(SQR(2))))
LET x = (a - 1)/(a + 1)
LET d = 1
LET s = x
LET a = x*x
FOR i = 1 TO k
   LET d = d + 2
   LET x = x*a
   LET s = s + x/d
NEXT i
LET s = 32*s
LET r = TIME - t
LET r$ = STR$(s)                           
PRINT
PRINT "0";
PRINT r$(0:1);
FOR i = 2 TO nd + 1
   PRINT r$(i:i);
   IF MOD((i - 1),10) = 0 THEN PRINT " ";
   IF MOD((i - 1),50) = 0 THEN 
      PRINT
      PRINT "  ";
   END IF
NEXT i
IF MOD (i - 2,50) <> 0  OR nd = 0 THEN PRINT
PRINT 
PRINT "Runtime: ";
PRINT  USING "0.##": r;
PRINT " seconds"
END

--------------------------------------------------- 

Log(2) =

0.6931471805 5994530941 7232121458 1765680755 0013436025 
  5254120680 0094933936 2196969471 5605863326 9964186875 
  4200148102 0570685733 6855202357 5813055703 2670751635 
  0759619307 2757082837 1435190307 0386238916 7347112335 
  0115364497 9552391204 7517268157 4932065155 5247341395 
  2588295045 3007095326 3666426541 0423915781 4952043740 
  4303855008 0194417064 1671518644 7128399681 7178454695 
  7026271631 0645461502 5720740248 1637773389 6385506952 
  6066834113 7273873722 9289564935 4702576265 2098859693 
  2019650585 5476470330 6793654432 5476327449 5125040606 
  9438147104 6899465062 2016772042 4524529612 6879465461 
  9316517468 1392672504 1038025462 5965686914 4192871608 
  2938031727 1436778265 4877566485 0856740776 4845146443 
  9940461422 6031930967 3540257444 6070308096 0850474866 
  3852313818 1676751438 6674766478 9088143714 1985494231 
  5199735488 0375165861 2753529166 1000710535 5824987941 
  4729509293 1138971559 9820565439 2871700072 1808576102 
  5236889213 2449713893 2037843935 3088774825 9701715591 
  0708823683 6275898425 8918535302 4363421436 7061189236 
  7891923723 1467232172 0534016492 5687274778 2344535343 
  
Runtime: 0.03 seconds

Thank you, Pauli and Mike, for your contributions!

Gerson.

Edited to fix a typo in line 7 of Mike’s first RPN program (RCL X, not RCL Y).
Also, 1 E in line 8 has been replaced with 1, as in his second program. That
is one byte shorter and doesn’t cause Free42 to show an “Internal Error”
message upon execution, when being pasted to it. Somehow Free42 3.1.10
doesn’t accept the constant 1E or 1 E, unless these are entered manually.