Creating digits of pi

[Gerson W. Barbosa]

(02-15-2018 11:57 PM)Gerson W. Barbosa Wrote:  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.

Sorry to post again in this relatively old thread, but I think it still fits here.

Yesterday, while googling for Madhava's correction terms (which are simply the first convergents of the continued fraction above, itself not unbeknownst to Hindu mathematicians of that time, according to some sources), I found a 5-hundred-year old infinite series for \(\pi\), by Nilakantha:

\(\frac{\pi}{8} = \frac{1}{1\cdot 3} + \frac{1}{5\cdot 7 } + \frac{1}{9\cdot 11} +\frac{1}{13\cdot 15}+\cdots =\sum_{n=1}^{\infty}\frac{1}{\left ( 4n-3 \right )\left ( 4n-1 \right )}\)

( Pride of India: A Glimpse Into India's Scientific Heritage, page 53 )

It's a small improvement of Madhava's original series (there are other two in the reference above). It can overwhelmingly be accelerated by a continuous fraction as well:

\(\frac{\pi}{8} \approx \frac{1}{1\cdot 3} + \frac{1}{5\cdot 7 } + \frac{1}{9\cdot 11} +\frac{1}{13\cdot15}+\cdots+\frac{1}{\left ( 4n-3 \right )\left ( 4n-1 \right )} +\frac{1}{16n+\frac{1}{n+\frac{4}{16n+\frac{9}{n+\frac{16}{\ddots +\frac{\ddots}{\ddots+\frac{n^{2}}{n\left [ 16-15\left ( n\bmod 2\right ) \right ]}}}}}}}\)


Being a non-alternating series, it is somewhat more simple to program, giving slightly more than two digits per iteration, when only one term of the series and one term of the continued fraction are processed at a time, starting from the respective last terms:


OPTION ARITHMETIC DECIMAL_HIGH       
INPUT  PROMPT "Number of decimal places: ":f                        
LET t = TIME                          
LET k = CEIL((f + 1)*12/25)
LET s = 0
LET a = 16*k*(k - 1) + 3  
LET b = 32*k
LET c = 0
LET n = k*k
LET m = 2*k - 1
LET d = k*(16 - 15*MOD(k,2))
LET e = 15*k*(2*MOD(k,2) - 1)
FOR i = 1 TO k
   LET s = s + 1/a
   LET b = b - 32
   LET a = a - b
   LET c = n/(c + d)
   LET n = n - m
   LET m = m - 2
   LET d = d + e
   LET e = -e
NEXT i 
LET c = 1/(c + d)
LET r = TIME - t                            
PRINT
LET p$ = STR$(8*(s + c))
PRINT p$(1:2);
FOR i = 3 TO f + 2
   PRINT p$(i:i);
   IF MOD((i - 2),10) = 0 THEN PRINT " ";
   IF MOD((i - 2),100) = 0 THEN 
      PRINT
      PRINT "  ";
   END IF
NEXT i
IF MOD (i - 3,100) <> 0  OR f = 0 THEN PRINT
PRINT
PRINT "Runtime: ";
PRINT  USING "0.##": r;
PRINT " seconds"
END



Number of decimal places: 999

3.1415926535 8979323846 2643383279 5028841971 6939937510 5820974944 5923078164 0628620899 8628034825 3421170679 
  8214808651 3282306647 0938446095 5058223172 5359408128 4811174502 8410270193 8521105559 6446229489 5493038196 
  4428810975 6659334461 2847564823 3786783165 2712019091 4564856692 3460348610 4543266482 1339360726 0249141273 
  7245870066 0631558817 4881520920 9628292540 9171536436 7892590360 0113305305 4882046652 1384146951 9415116094 
  3305727036 5759591953 0921861173 8193261179 3105118548 0744623799 6274956735 1885752724 8912279381 8301194912 
  9833673362 4406566430 8602139494 6395224737 1907021798 6094370277 0539217176 2931767523 8467481846 7669405132 
  0005681271 4526356082 7785771342 7577896091 7363717872 1468440901 2249534301 4654958537 1050792279 6892589235 
  4201995611 2129021960 8640344181 5981362977 4771309960 5187072113 4999999837 2978049951 0597317328 1609631859 
  5024459455 3469083026 4252230825 3344685035 2619311881 7101000313 7838752886 5875332083 8142061717 7669147303 
  5982534904 2875546873 1159562863 8823537875 9375195778 1857780532 1712268066 1300192787 6611195909 216420198

Runtime: 0.30 seconds

The 12/25 factor to determine the required number of terms may not be exact. I expected it to be 1/2 from the first tests, but I cannot confirm it to be correct the thousand-digit evaluation which is possible in Decimal BASIC.

Edited to remove Text Size large option from one of the formulas (which was stopping the page from opening is some Apple devices with iOS 11)
Edited to replace mid$ and left$ with their ANSI equivalents, per Tom L's suggestion below. Also, output format has been changed to display digits in groups of ten.


06-01-2018, 12:27 AM

PS:

(HP 50g) RPL versions:

1) Built-in LongInt version

%%HP: T(3)A(R)F(.);
\<< R\->I PUSH RAD -105 CF -3 CF DUP 1 + 12 * 25 / CEIL DUPDUP 1 - * 16 * 3 + OVER 32 * ROT DUPDUP SQ UNROT DUP + 1 - OVER DUP 2 MOD 15 * NEG 16 + * ROT DUPDUP 2 MOD 2 * 1 - * 15 * 0 8 ROLLD 7 ROLL ROT 0 UNROT 8 ROLL 1 ROT
START UNROT DUP INV ROT + ROT 32 - ROT OVER - 7 PICK 9 ROLL 7 PICK + / 8 ROLLD 7 ROLL 7 PICK - 7 ROLLD 6 ROLL 2 - 6 ROLLD 5 ROLL 5 PICK + 5 ROLLD 4 ROLL NEG 4 ROLLD SWAP
NEXT 8 ROLL 6 ROLL + INV 4 ROLL + 6 ROLLD 5 DROPN 8 * EXPAND FXND PICK3 ALOG OVER - PICK3 * SWAP IQUOT + \->STR DUP HEAD -51 FC? { "." } { "," } IFTE + SWAP TAIL + 1 ROT 2 + SUB POP
\>>

It will become extremely slow for increasing number of digits:

100 digits:    54.77 seconds
200 digits:   437.91 seconds
500 digits: 20583.40 seconds

2) LongFloat version

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

 100 digits:   24.98 seconds
 200 digits:   83.11 seconds
 500 digits:  731.91 seconds
1000 digits: 5111.75 seconds

2000 digits:  851.03 seconds (on the emulator, @ 2.6 GHz)

3.
14159265358979323846264338327950288419716939937510
58209749445923078164062862089986280348253421170679
82148086513282306647093844609550582231725359408128
48111745028410270193852110555964462294895493038196
44288109756659334461284756482337867831652712019091
45648566923460348610454326648213393607260249141273
72458700660631558817488152092096282925409171536436
78925903600113305305488204665213841469519415116094
33057270365759591953092186117381932611793105118548
07446237996274956735188575272489122793818301194912
98336733624406566430860213949463952247371907021798
60943702770539217176293176752384674818467669405132
00056812714526356082778577134275778960917363717872
14684409012249534301465495853710507922796892589235
42019956112129021960864034418159813629774771309960
51870721134999999837297804995105973173281609631859
50244594553469083026425223082533446850352619311881
71010003137838752886587533208381420617177669147303
59825349042875546873115956286388235378759375195778
18577805321712268066130019278766111959092164201989
38095257201065485863278865936153381827968230301952
03530185296899577362259941389124972177528347913151
55748572424541506959508295331168617278558890750983
81754637464939319255060400927701671139009848824012
85836160356370766010471018194295559619894676783744
94482553797747268471040475346462080466842590694912
93313677028989152104752162056966024058038150193511
25338243003558764024749647326391419927260426992279
67823547816360093417216412199245863150302861829745
55706749838505494588586926995690927210797509302955
32116534498720275596023648066549911988183479775356
63698074265425278625518184175746728909777727938000
81647060016145249192173217214772350141441973568548
16136115735255213347574184946843852332390739414333
45477624168625189835694855620992192221842725502542
56887671790494601653466804988627232791786085784383
82796797668145410095388378636095068006422512520511
73929848960841284886269456042419652850222106611863
06744278622039194945047123713786960956364371917287
46776465757396241389086583264599581339047802759010