Creating digits of pi

[Gerson W. Barbosa]

(02-14-2018 11:57 PM)Gerson W. Barbosa Wrote:  By using this method, a thousand digits can be computed in less than four hours (3h51m09s) on a 49G.

This version is slightly faster, 3h44m25s on the 49G (1h58m17s on the 50g, 2 minutes on the emulator).

Code:


%%HP: T(3)A(D)F(.);
\<< 1. + DUP 'DIGITS' STO 4. + 3. / IP DUPDUP + R\<-\->F OVER SQ 4. * 
OVER R\<-\->F DUP + DUPDUP + R\<-\->F SWAP 1. - 0. R\<-\->F DUP 6. ROLLD 
1. 8. ROLL 
  START OVER R\<-\->F FINV FSUB OVER 2. - R\<-\->F FINV FADD 6. ROLLD ROT 
DUP2 SWAP - R\<-\->F 5. PICK PICK3 R\<-\->F 6. PICK 9. ROLL FADD FDIV 
FADD FDIV 6. ROLLD OVER 1. - DUP + - UNROT 4. - 5. ROLL 
  NEXT ROT 6. ROLL FADD FINV FADD 4. R\<-\->F FMULT 4. ROLLD 3. DROPN 
ZZ\<-\->F DROP \->STR DUP HEAD -51. FC? { "." } { "," } IFTE + SWAP TAIL + 
\>>



d = 1000

3.1415926535897932384626433832795028841971693993751058209749445923078164062
862089986280348253421170679821480865132823066470938446095505822317253594081
284811174502841027019385211055596446229489549303819644288109756659334461284
756482337867831652712019091456485669234603486104543266482133936072602491412
737245870066063155881748815209209628292540917153643678925903600113305305488
204665213841469519415116094330572703657595919530921861173819326117931051185
480744623799627495673518857527248912279381830119491298336733624406566430860
213949463952247371907021798609437027705392171762931767523846748184676694051
320005681271452635608277857713427577896091736371787214684409012249534301465
495853710507922796892589235420199561121290219608640344181598136297747713099
605187072113499999983729780499510597317328160963185950244594553469083026425
223082533446850352619311881710100031378387528865875332083814206171776691473
035982534904287554687311595628638823537875937519577818577805321712268066130
019278766111959092164201989

  :s:13479.2559 (HP 49g)

   :s:7097.2896 (HP 50g)

    :s:120.5122 (Emu50g @ 2.6 GHz)

Formula:

\(\frac{\pi }{4}\approx 1-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\cdots +\frac{1}{2n-3}-\frac{1}{2n-1}+\frac{1}{4n+\frac{1^{2}}{n+\frac{2^{2}}{4n+\frac{3^{2}}{n+\frac{4 ^{2}}{...+\frac{...}{...+\frac{n^{2}}{4n}}}}}}}\)


This gives \(\frac{3n}{2}\) correct digits (for even n).

Algorithm:

 d = 1000                                       ; desired number of digits
 n = int(d/3) + 1                               ; required number of iterations
 a = 2*n                                        ; number of terms in the series and in the continued fractions on (n, in the above formula)
 b = 8*n                                        ; part of the denominators in the cont'ed fraction, alternating with the previous constant (4n and n, in the formula)
 c = 4*n*n                                      ; rightmost squared integer in the continued fraction (n² in the formula)
 d = 4*n - 1                                    ; rightmost denominator in the series (n-1 in the formula)  
 s1 = 0                                         ; initialization of the first partial sum (terms of series)
 s2 = 0                                         ; initialization of the second partial sum (terms of the continued fraction)
 for i = 1 to n                                 ; start loop
    s1 = s1 - 1/d + 1/(d - 2)                   ; both the terms of the series and the terms of the continued fraction are    
    s2 = (c - d)/(a + c/(b + s2))               ; processed two at a time, starting from the respective rightmost terms
    c = c - 2*(d - 1)                           ; next squared integers in the continued fraction are computed by simple subtraction,     
    d = d - 4                                   ; like the the next denominators in the series
 next i                                         ; end loop
 s2 = 1/(b + s2)                                ; add leftmost term of the continued fraction to the respective partial sum
 pi = 4*(s1 + s2)                               ; add both partial sums together and multiply the result by 4 to get pi
 display pi                                     ; display pi to d decimal places


This gives \(3\) digits per iteration.

Edited to fix a couple of typos

Edited again to fix another typo per Mike’s message below