Leibniz formula for π on HP 35s

[lipoff]

I got inspired by this vintage computer race to use the slow-converging alternating Leibnitz series method to approximate π.

It's a few years old, but in the video an HP 9825B from 1980 takes 60 seconds to compute π to four decimal places, while an Android Fairphone 2 from 2015 takes 6 seconds to make the same computation. The HP was programmed in HPL; the Android in Python. Obviously the 2.26 GHz droid could run a lot faster, but it's running an interpreted language with lots of overhead.

How well can the underpowered 33 kHz processor in the HP 35s perform in RPN? My first attempt took 4 minutes and 45 seconds, but this program, my fifth revision, is 28 instructions long and takes 2 minutes and 40 seconds (160 seconds) to calculate 4 digits of π on the HP 35s:

Code:


L001    LBL L
L002    0
L003    STO P
L004    1
L005    STO L
L006    0.84501
L007    STO Z
L008    RCL L
L009    1/x
L010    RCL L
L011    2
L012    +
L013    1/x
L014    -
L015    STO+ P
L016    4
L017    STO+ L
L018    ISG Z
L019    GTO L008
L020    RCL P
L021    0.78525
L022    x>y?
L023    GTO L006
L024    4
L025    RCLx P
L026    2
L027    RCL+ L
L028    RTN

Controlling for clockspeed, that's 68x more efficient than the HP 9825B and 2750x more efficient than the smartphone! Of course, I'm cheating a little by running four loops of precisely 845 iterations, but even with a less arbitrary/optimized value of 1000, it's almost as fast.

[PedroLeiva]

Very interesting, I have HP 35s. Could you please include instructions to run this, and an example to beging, Pedro
Just pressing XEQ L ENTER I got:
3.141001637
3.387.0000000

[Gerson W. Barbosa]

(05-06-2020 02:45 PM)lipoff Wrote:  I got inspired by this vintage computer race to use the slow-converging alternating Leibnitz series method to approximate π.

It’s now known as Madhava-Leibniz Series, as in the video (1:09). Madhava himself would never compute π this way, even if he had a fast computer. Anyway, that’s an interesting experiment.
Madhava would compute just a few terms then he would apply his best correction term to get more digits. Those correction terms (he has provided three of them) can be directly obtained from a continued fraction that was probably known to him. See Madhava of Sangamagrama – The value of π.

Here is an HP-75C program that computes π using the series and the continued fraction. (Commented listing here).

The program asks for the number of digits, but at this early stage an input of 9 will give the familiar 12-digit value:

1 INPUT D
2 N=IP(D/3)+1
3 A=2*N
4 B=8*N
5 C=4*N*N
6 D=4*N-1
7 S1=0
8 S2=0
9 FOR I=1 TO N
10 S1=S1-1/D+1/(D-2)
11 S2=(C-D)/(A+C/(B+S2))
12 C=C-2*(D-1)
13 D=D-4
14 NEXT I
15 S2=1/(B+S2)
16 DISP 4*(S1+S2)

RUN
?9 RTN
-> 3.14159265359 (4 iterations, 0.304 seconds)

As a side note, checking the catalog I noticed two unintentional coincidences:

CAT
-> PI B 314 12:56 06\05\20

314 bytes (~100π), file created at 12:56 hour (~4π) :-)

[lipoff]

(05-06-2020 05:09 PM)PedroLeiva Wrote:  Very interesting, I have HP 35s. Could you please include instructions to run this, and an example to beging, Pedro
Just pressing XEQ L ENTER I got:
3.141001637
3.387.0000000

Yes, you ran it just right. It computed π to four decimal places, and 1/3387 was the last term in the series.

It's a terrible way to compute π from a computational efficiency standpoint although remarkable that such an alternating series sums to π/4 and it was fun to implement.

If I understand ISG and DSE correctly on the HP 35s the maximum number of steps for an ISG loop is 10^3 (Z = 0.99901) but for a DSE loop is 10^7 (Z = 9999999.00001). In this case it doesn't matter; 844 steps is enough.

[lipoff]

If it's helpful, here's the program with comments, with a very slight update to make it "more optimum" as it now runs in 2 minutes and 34 seconds:

Code:

L001    LBL L                   Name the program
L002    0            
L003    STO P                   Set our accumulator for the series to 0
L004    1
L005    STO L                   Start the series denominator at 1
L006    1.84401            
L007    STO Z                   Start a loop counter at 1, with a max of 844 and increment of 1
L008    RCL L                   This is the start of our calculation loop
L009    1/x                     Take 1/L
L010    STO+ P                  Add this quantity to P
L011    RCL L 
L012    2
L013    +
L014    1/x                     Take 1/(L+2)
L015    STO- P                  Subtract this quantity from P
L016    4
L017    STO+ L                  Add 4 to our denominator
L018    ISG Z                   Increment the loop counter and check if it exceeds the max value
L019    GTO L008                If not, loop back to L008
L020    RCL P            
L021    0.78525            
L022    x>y?                    Check if our estimate, P < 3.141/4 = 0.78525 because the series converges to pi/4
L023    GTO L006                If so, loop back to L006
L024    4
L025    RCLx P                  If not, display 4 times P    
L026    2
L027    RCL+ L                  And display the final denominator in the series
L028    RTN                     End program