Creating digits of pi

[Gerson W. Barbosa]

(02-10-2018 12:44 PM)Mike (Stgt) Wrote:  Could someone of the interested (Decimal BASIC admirer, etc.) please confirm this Ramanujan–Sato of 1993 works? I fail with it, get Pi = 4.1748739, what is a bit too far off.
I assume a typo in the Wiki.

Ciao.....Mike

Off by a factor of exactly 400/301:

4.1748739 * 301/400 = 3.1415926

Probably just a coincidence, but I would check this out with more digits.

[toml_12953]

(02-10-2018 12:44 PM)Mike (Stgt) Wrote:  Could someone of the interested (Decimal BASIC admirer, etc.) please confirm this Ramanujan–Sato of 1993 works? I fail with it, get Pi = 4.1748739, what is a bit too far off.
I assume a typo in the Wiki.

Ciao.....Mike

The quadratic convergence works:

Code:

OPTION ARITHMETIC DECIMAL_HIGH
LET a=SQR(2)
LET b=0
LET p=2+SQR(2)
LET n=1
DO UNTIL a=1
   LET an=a
   LET bn=b
   LET pn=p
   LET a=(SQR(an)+1/SQR(an))/2
   LET b=((1+bn)*SQR(an))/(an+bn)
   LET p=((1+a)*pn*b)/(1+b)
LOOP
PRINT p
END

[toml_12953]

(02-10-2018 04:54 PM)Mike (Stgt) Wrote:  

(02-10-2018 04:42 PM)toml_12953 Wrote:  ... Something is off there somewhere.

Yes, have a closer look to the indices.

Try this:

Code:

an = sqrt(2)
bn = 0
pn = 2 + an
do until p = pn
   a = an
   b = bn
   p = pn
   ar = sqrt(a)
   an = .5 * (ar + 1 / ar)
   bn = ar * (1 + b) / (a + b)
   pn = p * bn * (1 + an) / (1 + bn)
end
say p    /* or the last 30 digits if it's too long */

an stands for a(n+1), a stands for a(n), and so on.

Ciao.....Mike

I corrected my original code. It works now. Thanks!

[Valentin Albillo]

(02-10-2018 12:44 PM)Mike (Stgt) Wrote:  Could someone of the interested (Decimal BASIC admirer, etc.) please confirm this Ramanujan–Sato of 1993 works? I fail with it, get Pi = 4.1748739, what is a bit too far off.
I assume a typo in the Wiki.

Nope, the Wiki is correct. I used the constants featured there and this is the result I got:

list

10 word -520:point -160
20 A=63365028312971999585426220+28337702140800842046825600*sqrt(5)+384*sqrt
(5)*sqrt(10891728551171178200467436212395209160385656017+48709290865788102250773​
38534541688721351255040*sqrt(5))
30 B=7849910453496627210289749000+3510586678260932028965606400*sqrt(5)+2515
968*sqrt(3110)*sqrt(6260208323789001636993322654444020882161+2799650273060444296​
577206890718825190235*sqrt(5))
40 C=-214772995063512240-96049403338648032*sqrt(5)-1296*sqrt(5)*sqrt(109852
34579463550323713318473+4912746253692362754607395912*sqrt(5))
50 clr time : S=0:N=0
60 T=S+!(6*N)/!(3*N)/(!(N)^3)*(A+N*B)/C^(3*N):if T<>S then S=T:inc N:goto 60
70 print sqrt((-C)^3)/S:print time1000


run

3.14159265358979323846264338327950288419716939937510582097494459230781640628620
89986280348253421170679821480865132823066470938446095505822317253594081284811174​
50284102701938521105559644622948954930381964428810975665933446128475648233786783​
16527120190914564856692346034861045432664821339360726024914127372458700660631558​
81748815209209628292540917153643678925903600113305305488204665213841469519415116​
09433057270365759591953092186117381932611793105118548074462379962749567351885752​
72489122793818301194912983367336244065664308602139494639522473719070217986094370​
27705392171762931767523846748184676694051320005681271452635608277857713427577896​
09173637178721468440901224953430146549585371050792279689258923542019956112129021​
96086403441815981362977477130996051870721134999999837

which gives 770 correct decimals (all featured) in 15-16 iterations, 8 milliseconds.

Have a nice weekend.
V.
.

[emece67]

I've tested, using Python + the mpmath package, the nonic convergence method. It gives

• 1 correct digit at the 1st iteration
• 22 at the 2nd
• 217 at the 3rd
• 1984 at the 4th
• 17898 at the 5th
• 161123 at the 6th
• 1450163 at the 7th
  

Really awesome!

Regards.

[TASP]

Yeah, not helpful at all, but for some random digit of Pi in base 2, just guess.

50/50 chance of being right is good enough for me and I don't need to put the batteries back in my HP41 to do it.

[Gerson W. Barbosa]

(02-09-2018 07:27 PM)Gerson W. Barbosa Wrote:  

Code:


OPTION ARITHMETIC DECIMAL_HIGH ! 1000 digits precision
LET t1=TIME
LET d=1000                     ! 1000 decimal places
LET n=INT(D/3) + 1
LET a=2*n
LET b=8*n  
LET s0=0
LET s1=0
FOR i=n TO 1 STEP -1
   LET s0=s0 + 2/(16*i*(i - 1) + 3)
   LET t=4*i
   LET s1=((b + s1)*((i - 1)*t + 1))/(a*(b + s1) + i*t)
NEXT i
LET s1=1/(b + s1)
LET p=4*(s0 + s1)
PRINT p
PRINT TIME - t1 
END

334 iterations (but 327 would have sufficed for 1000 decimal places), 342 ms (average of 5 runs). Still slower than Tom's version (243 ms here, Intel E4700 @ 2.6 GHz), but about 130 % faster when compared to my previous version above).

Now I get running times as lows as 0.2 seconds, but sometimes as high as 0.35 seconds.

Code:


OPTION ARITHMETIC DECIMAL_HIGH ! 1000 digits precision
LET t = TIME
LET d = 1000                   ! 1000 decimal places
LET n = INT(D/3+1/2)
LET a = 2*n
LET b = 8*n 
LET c = 4*n*n
LET d = 4*n - 1 
LET s0 = 0
LET s1 = 0
FOR i = 1 TO n
   LET s0 = s0 + 2/((d - 2)*d)
   LET s1 = (b + s1)*(c - d)/(a*(b + s1) + c)
   LET c = c - d - d + 2   
   LET d = d - 4
NEXT i
LET s1 = 1/(b + s1)
LET p=4*(s0 + s1)
PRINT p
PRINT
PRINT TIME - t; "seconds"  
END

3.141592653589793238462643383279502884197169399375105820974944592307816406286208​99862803482534211706798214808651328230664709384460955058223172535940812848111745​02841027019385211055596446229489549303819644288109756659334461284756482337867831​65271201909145648566923460348610454326648213393607260249141273724587006606315588​17488152092096282925409171536436789259036001133053054882046652138414695194151160​94330572703657595919530921861173819326117931051185480744623799627495673518857527​24891227938183011949129833673362440656643086021394946395224737190702179860943702​77053921717629317675238467481846766940513200056812714526356082778577134275778960​91736371787214684409012249534301465495853710507922796892589235420199561121290219​60864034418159813629774771309960518707211349999998372978049951059731732816096318​59502445945534690830264252230825334468503526193118817101000313783875288658753320​83814206171776691473035982534904287554687311595628638823537875937519577818577805​32171226806613001927876611195909216420199

.2 seconds


These are the more straightforward expressions for s0 and s1, but they're slower:

Code:


LET s0 = s0 - 1/d + 1/(d - 2)
LET s1 = (c - d)/(a + c/(b + s1))

[Gerson W. Barbosa]

(02-12-2018 01:38 AM)Mike (Stgt) Wrote:  BTW, 1000 decimals in 10 minutes on an HP-49, that gives everyone the chance to keep count of it.

Really fast (and compact), 5 seconds on the emulator!

For real fun, I would try a 2500-year-old algorithm (with a small adjustment for speed) on my 50g. No problem waiting three hours for the glorious thousand digits :-)

http://www.hpmuseum.org/cgi-sys/cgiwrap/...443#189221

There is also a 4-minute alternative (certainly more than 10 minutes on the 49G). Both programs require the LongFloat library, though.

[ttw]

Some methods for decimal digits but not quite as nice as the base 16 version.

http://numbers.computation.free.fr/Const...ldigit.pdf

https://arxiv.org/ftp/arxiv/papers/0912/0912.0303.pdf

Makes a good test problem for compilers and hardware.

[Massimo Gnerucci]

(02-14-2018 04:55 AM)Mike (Stgt) Wrote:  Just tried on an HP300s+

ln(640320)/sqrt(163)

and got as result Pi/3 -- very impressive, seems as if this calculator could compute Pi to the last digit.

Ciao.....Mike

There's no last digit in π... ;)

[Joe Horn]

(02-14-2018 09:52 AM)Mike (Stgt) Wrote:  BTW, it must be a good approximation so the 300s+ will take it as Pi, 355/113 is not good enough.

It's very good. ln(640320)/sqrt(163) differs from pi/3 by only 7.39889E-17 (rounded of course). The simplest fraction that's even closer to pi/3 is 165707065/158238591 (error is 5.46945E-17). It's amazing what an insomniac brain will occupy itself with at 3:06AM.

[Massimo Gnerucci]

(02-14-2018 09:52 AM)Mike (Stgt) Wrote:  

(02-14-2018 07:37 AM)Massimo Gnerucci Wrote:  There's no last digit in π...

This 'last digit of Pi' was only a reference to Star Trek
Kirk asks: "Aren't there some mathematical problems that simply can't be solved?"
And Spock 'fries the brains' of a rogue computer by telling it: "Compute to the last digit the value of Pi."



BTW, it must be a good approximation so the 300s+ will take it as Pi, 355/113 is not good enough.

/M.

Oh yes! Wolf in the fold.

[John Keith]

(02-12-2018 03:40 AM)Gerson W. Barbosa Wrote:  

(02-12-2018 01:38 AM)Mike (Stgt) Wrote:  BTW, 1000 decimals in 10 minutes on an HP-49, that gives everyone the chance to keep count of it.

Really fast (and compact), 5 seconds on the emulator!

For real fun, I would try a 2500-year-old algorithm (with a small adjustment for speed) on my 50g. No problem waiting three hours for the glorious thousand digits :-)

http://www.hpmuseum.org/cgi-sys/cgiwrap/...443#189221

There is also a 4-minute alternative (certainly more than 10 minutes on the 49G). Both programs require the LongFloat library, though.

I know it's "cheating", but LongFloat's built-in function FPI returns 1000 digits of PI in about 42 seconds on the HP50, about 1.1s on the emulator. :-)

[Gerson W. Barbosa]

(02-14-2018 07:24 PM)John Keith Wrote:  

(02-12-2018 03:40 AM)Gerson W. Barbosa Wrote:  For real fun, I would try a 2500-year-old algorithm (with a small adjustment for speed) on my 50g. No problem waiting three hours for the glorious thousand digits :-)

http://www.hpmuseum.org/cgi-sys/cgiwrap/...443#189221

There is also a 4-minute alternative (certainly more than 10 minutes on the 49G). Both programs require the LongFloat library, though.

I know it's "cheating", but LongFloat's built-in function FPI returns 1000 digits of PI in about 42 seconds on the HP50, about 1.1s on the emulator. :-)

Less than 74 seconds on my HP 49G. But this is not about getting it fast, it's about testing a method that was available in India at least sinced the late 14th century. I believe the continued fraction which works as a correction to the very slowly convergent Madhava-Leibniz series (a.k.a. Gregory series), was known to Madhava. 20 terms of the series and 20 terms of the continued fraction would have yielded 30 correct digits, but the latter might be not practical for manual calculation, or they just thought it was not worth the effort.

By using this method, a thousand digits can be computed in less than four hours (3h51m09s) on a 49G. I don't think they would have found that too slow back then :-)

1000

Code:


%%HP: T(3)A(D)F(.);
\<< 1. + DUP 'DIGITS' STO 4. + 3. / IP R\<-\->F DUPDUP FADD OVER FSQ 4. R\<-\->F 
FMULT OVER DUP FADD DUPDUP FADD SWAP 1. R\<-\->F FSUB 0. R\<-\->F DUP 1. 8. ROLL 
R\<-\->F
  START UNROT OVER FINV FSUB OVER 2. R\<-\->F FSUB FINV FADD 6. ROLLD 4. ROLLD 
OVER FADD PICK3 SWAP FDIV 4. ROLLD UNROT DUP2 SWAP FSUB 5. ROLL 6. PICK FADD FDIV
 UNROT OVER 1. R\<-\->F FSUB DUP FADD FSUB UNROT 4. R\<-\->F FSUB UNROT 5. ROLLD 
UNROT 6. ROLL 6. ROLL
  NEXT 4. PICK FADD FINV FADD 4. R\<-\->F FMULT 5. ROLLD 4. DROPN
\>>

EVAL

-->


 31415926535897932384626433832795028841971693993751058209749445923078164062862
089986280348253421170679821480865132823066470938446095505822317253594081284811
174502841027019385211055596446229489549303819644288109756659334461284756482337
867831652712019091456485669234603486104543266482133936072602491412737245870066
063155881748815209209628292540917153643678925903600113305305488204665213841469
519415116094330572703657595919530921861173819326117931051185480744623799627495
673518857527248912279381830119491298336733624406566430860213949463952247371907
021798609437027705392171762931767523846748184676694051320005681271452635608277
857713427577896091736371787214684409012249534301465495853710507922796892589235
420199561121290219608640344181598136297747713099605187072113499999983729780499
510597317328160963185950244594553469083026425223082533446850352619311881710100
031378387528865875332083814206171776691473035982534904287554687311595628638823
537875937519577818577805321712268066130019278766111959092164201989.E-1000

s: 13868.7266

That's about the time it takes the HP 50g to get the equivalent algebraic expression. Expanding it would take much, much longer.

1000

Code:


« PUSH RAD -105 CF -3 CF 4 + 3 IDIV2 DROP DUPDUP + OVER SQ 4 * OVER DUP + DUPDUP + SWAP 1 - 0 DUP 1
 8 ROLL
  START UNROT OVER INV - OVER 2 - INV + 6 ROLLD 4 ROLLD OVER + PICK3 SWAP / 4 ROLLD UNROT DUP2
 SWAP - 5 ROLL 6 PICK + / UNROT OVER 1 - DUP + - UNROT 4 - UNROT 5 ROLLD UNROT 6 ROLL 6 ROLL
  NEXT 4
 PICK + INV + 4 * 5 ROLLD 4 DROPN POP
»

EVAL

-->

Code:


'(-(1/1335)+1/1333-1/1331+1/1329-1/1327+1/1325-1/1323+1/1321-1/1319+1/1317-1/1315+1/1313-1/1311+1/1309-
1/1307+1/1305-1/1303+1/1301-1/1299+1/1297-1/1295+1/1293-1/1291+1/1289-1/1287+1/1285-1/1283+1/1281-
1/1279+1/1277-1/1275+1/1273-1/1271+1/1269-1/1267+1/1265-1/1263+1/1261-1/1259+1/1257-1/1255+1/1253-
1/1251+1/1249-1/1247+1/1245-1/1243+1/1241-1/1239+1/1237-1/1235+1/1233-1/1231+1/1229-1/1227+1/1225-
1/1223+1/1221-1/1219+1/1217-1/1215+1/1213-1/1211+1/1209-1/1207+1/1205-1/1203+1/1201-1/1199+1/1197-
1/1195+1/1193-1/1191+1/1189-1/1187+1/1185-1/1183+1/1181-1/1179+1/1177-1/1175+1/1173-1/1171+1/1169-
1/1167+1/1165-1/1163+1/1161-1/1159+1/1157-1/1155+1/1153-1/1151+1/1149-1/1147+1/1145-1/1143+1/1141-
1/1139+1/1137-1/1135+1/1133-1/1131+1/1129-1/1127+1/1125-1/1123+1/1121-1/1119+1/1117-1/1115+1/1113-
1/1111+1/1109-1/1107+1/1105-1/1103+1/1101-1/1099+1/1097-1/1095+1/1093-1/1091+1/1089-1/1087+1/1085-
1/1083+1/1081-1/1079+1/1077-1/1075+1/1073-1/1071+1/1069-1/1067+1/1065-1/1063+1/1061-1/1059+1/1057-
1/1055+1/1053-1/1051+1/1049-1/1047+1/1045-1/1043+1/1041-1/1039+1/1037-1/1035+1/1033-1/1031+1/1029-
1/1027+1/1025-1/1023+1/1021-1/1019+1/1017-1/1015+1/1013-1/1011+1/1009-1/1007+1/1005-1/1003+1/1001-
1/999+1/997-1/995+1/993-1/991+1/989-1/987+1/985-1/983+1/981-1/979+1/977-1/975+1/973-1/971+1/969-
1/967+1/965-1/963+1/961-1/959+1/957-1/955+1/953-1/951+1/949-1/947+1/945-1/943+1/941-1/939+1/937-
1/935+1/933-1/931+1/929-1/927+1/925-1/923+1/921-1/919+1/917-1/915+1/913-1/911+1/909-1/907+1/905-
1/903+1/901-1/899+1/897-1/895+1/893-1/891+1/889-1/887+1/885-1/883+1/881-1/879+1/877-1/875+1/873-
1/871+1/869-1/867+1/865-1/863+1/861-1/859+1/857-1/855+1/853-1/851+1/849-1/847+1/845-1/843+1/841-
1/839+1/837-1/835+1/833-1/831+1/829-1/827+1/825-1/823+1/821-1/819+1/817-1/815+1/813-1/811+1/809-
1/807+1/805-1/803+1/801-1/799+1/797-1/795+1/793-1/791+1/789-1/787+1/785-1/783+1/781-1/779+1/777-
1/775+1/773-1/771+1/769-1/767+1/765-1/763+1/761-1/759+1/757-1/755+1/753-1/751+1/749-1/747+1/745-
1/743+1/741-1/739+1/737-1/735+1/733-1/731+1/729-1/727+1/725-1/723+1/721-1/719+1/717-1/715+1/713-
1/711+1/709-1/707+1/705-1/703+1/701-1/699+1/697-1/695+1/693-1/691+1/689-1/687+1/685-1/683+1/681-
1/679+1/677-1/675+1/673-1/671+1/669-1/667+1/665-1/663+1/661-1/659+1/657-1/655+1/653-1/651+1/649-
1/647+1/645-1/643+1/641-1/639+1/637-1/635+1/633-1/631+1/629-1/627+1/625-1/623+1/621-1/619+1/617-
1/615+1/613-1/611+1/609-1/607+1/605-1/603+1/601-1/599+1/597-1/595+1/593-1/591+1/589-1/587+1/585-
1/583+1/581-1/579+1/577-1/575+1/573-1/571+1/569-1/567+1/565-1/563+1/561-1/559+1/557-1/555+1/553-
1/551+1/549-1/547+1/545-1/543+1/541-1/539+1/537-1/535+1/533-1/531+1/529-1/527+1/525-1/523+1/521-
1/519+1/517-1/515+1/513-1/511+1/509-1/507+1/505-1/503+1/501-1/499+1/497-1/495+1/493-1/491+1/489-
1/487+1/485-1/483+1/481-1/479+1/477-1/475+1/473-1/471+1/469-1/467+1/465-1/463+1/461-1/459+1/457-
1/455+1/453-1/451+1/449-1/447+1/445-1/443+1/441-1/439+1/437-1/435+1/433-1/431+1/429-1/427+1/425-
1/423+1/421-1/419+1/417-1/415+1/413-1/411+1/409-1/407+1/405-1/403+1/401-1/399+1/397-1/395+1/393-
1/391+1/389-1/387+1/385-1/383+1/381-1/379+1/377-1/375+1/373-1/371+1/369-1/367+1/365-1/363+1/361-
1/359+1/357-1/355+1/353-1/351+1/349-1/347+1/345-1/343+1/341-1/339+1/337-1/335+1/333-1/331+1/329-
1/327+1/325-1/323+1/321-1/319+1/317-1/315+1/313-1/311+1/309-1/307+1/305-1/303+1/301-1/299+1/297-
1/295+1/293-1/291+1/289-1/287+1/285-1/283+1/281-1/279+1/277-1/275+1/273-1/271+1/269-1/267+1/265-
1/263+1/261-1/259+1/257-1/255+1/253-1/251+1/249-1/247+1/245-1/243+1/241-1/239+1/237-1/235+1/233-
1/231+1/229-1/227+1/225-1/223+1/221-1/219+1/217-1/215+1/213-1/211+1/209-1/207+1/205-1/203+1/201-
1/199+1/197-1/195+1/193-1/191+1/189-1/187+1/185-1/183+1/181-1/179+1/177-1/175+1/173-1/171+1/169-
1/167+1/165-1/163+1/161-1/159+1/157-1/155+1/153-1/151+1/149-1/147+1/145-1/143+1/141-1/139+1/137-
1/135+1/133-1/131+1/129-1/127+1/125-1/123+1/121-1/119+1/117-1/115+1/113-1/111+1/109-1/107+1/105-
1/103+1/101-1/99+1/97-1/95+1/93-1/91+1/89-1/87+1/85-1/83+1/81-1/79+1/77-1/75+1/73-1/71+1/69-1/67+1/65-
1/63+1/61-1/59+1/57-1/55+1/53-1/51+1/49-1/47+1/45-1/43+1/41-1/39+1/37-1/35+1/33-1/31+1/29-1/27+1/25-
1/23+1/21-1/19+1/17-1/15+1/13-1/11+1/9-1/7+1/5-1/3+1+1/(1/(4/(9/(16/(25/(36/(49/(64/(81/(100/(121/(144/
(169/(196/(225/(256/(289/(324/(361/(400/(441/(484/(529/(576/(625/(676/(729/(784/(841/(900/(961/(1024/
(1089/(1156/(1225/(1296/(1369/(1444/(1521/(1600/(1681/(1764/(1849/(1936/(2025/(2116/(2209/(2304/(2401/
(2500/(2601/(2704/(2809/(2916/(3025/(3136/(3249/(3364/(3481/(3600/(3721/(3844/(3969/(4096/(4225/(4356/
(4489/(4624/(4761/(4900/(5041/(5184/(5329/(5476/(5625/(5776/(5929/(6084/(6241/(6400/(6561/(6724/(6889/
(7056/(7225/(7396/(7569/(7744/(7921/(8100/(8281/(8464/(8649/(8836/(9025/(9216/(9409/(9604/(9801/(10000/
(10201/(10404/(10609/(10816/(11025/(11236/(11449/(11664/(11881/(12100/(12321/(12544/(12769/(12996/
(13225/(13456/(13689/(13924/(14161/(14400/(14641/(14884/(15129/(15376/(15625/(15876/(16129/(16384/
(16641/(16900/(17161/(17424/(17689/(17956/(18225/(18496/(18769/(19044/(19321/(19600/(19881/(20164/
(20449/(20736/(21025/(21316/(21609/(21904/(22201/(22500/(22801/(23104/(23409/(23716/(24025/(24336/
(24649/(24964/(25281/(25600/(25921/(26244/(26569/(26896/(27225/(27556/(27889/(28224/(28561/(28900/
(29241/(29584/(29929/(30276/(30625/(30976/(31329/(31684/(32041/(32400/(32761/(33124/(33489/(33856/
(34225/(34596/(34969/(35344/(35721/(36100/(36481/(36864/(37249/(37636/(38025/(38416/(38809/(39204/
(39601/(40000/(40401/(40804/(41209/(41616/(42025/(42436/(42849/(43264/(43681/(44100/(44521/(44944/
(45369/(45796/(46225/(46656/(47089/(47524/(47961/(48400/(48841/(49284/(49729/(50176/(50625/(51076/
(51529/(51984/(52441/(52900/(53361/(53824/(54289/(54756/(55225/(55696/(56169/(56644/(57121/(57600/
(58081/(58564/(59049/(59536/(60025/(60516/(61009/(61504/(62001/(62500/(63001/(63504/(64009/(64516/
(65025/(65536/(66049/(66564/(67081/(67600/(68121/(68644/(69169/(69696/(70225/(70756/(71289/(71824/
(72361/(72900/(73441/(73984/(74529/(75076/(75625/(76176/(76729/(77284/(77841/(78400/(78961/(79524/
(80089/(80656/(81225/(81796/(82369/(82944/(83521/(84100/(84681/(85264/(85849/(86436/(87025/(87616/
(88209/(88804/(89401/(90000/(90601/(91204/(91809/(92416/(93025/(93636/(94249/(94864/(95481/(96100/
(96721/(97344/(97969/(98596/(99225/(99856/(100489/(101124/(101761/(102400/(103041/(103684/(104329/
(104976/(105625/(106276/(106929/(107584/(108241/(108900/(109561/(110224/(110889/(111556/(112225/
(112896/(113569/(114244/(114921/(115600/(116281/(116964/(117649/(118336/(119025/(119716/(120409/
(121104/(121801/(122500/(123201/(123904/(124609/(125316/(126025/(126736/(127449/(128164/(128881/
(129600/(130321/(131044/(131769/(132496/(133225/(133956/(134689/(135424/(136161/(136900/(137641/
(138384/(139129/(139876/(140625/(141376/(142129/(142884/(143641/(144400/(145161/(145924/(146689/
(147456/(148225/(148996/(149769/(150544/(151321/(152100/(152881/(153664/(154449/(155236/(156025/
(156816/(157609/(158404/(159201/(160000/(160801/(161604/(162409/(163216/(164025/(164836/(165649/
(166464/(167281/(168100/(168921/(169744/(170569/(171396/(172225/(173056/(173889/(174724/(175561/
(176400/(177241/(178084/(178929/(179776/(180625/(181476/(182329/(183184/(184041/(184900/(185761/
(186624/(187489/(188356/(189225/(190096/(190969/(191844/(192721/(193600/(194481/(195364/(196249/
(197136/(198025/(198916/(199809/(200704/(201601/(202500/(203401/(204304/(205209/(206116/(207025/
(207936/(208849/(209764/(210681/(211600/(212521/(213444/(214369/(215296/(216225/(217156/(218089/
(219024/(219961/(220900/(221841/(222784/(223729/(224676/(225625/(226576/(227529/(228484/(229441/
(230400/(231361/(232324/(233289/(234256/(235225/(236196/(237169/(238144/(239121/(240100/(241081/
(242064/(243049/(244036/(245025/(246016/(247009/(248004/(249001/(250000/(251001/(252004/(253009/
(254016/(255025/(256036/(257049/(258064/(259081/(260100/(261121/(262144/(263169/(264196/(265225/
(266256/(267289/(268324/(269361/(270400/(271441/(272484/(273529/(274576/(275625/(276676/(277729/
(278784/(279841/(280900/(281961/(283024/(284089/(285156/(286225/(287296/(288369/(289444/(290521/
(291600/(292681/(293764/(294849/(295936/(297025/(298116/(299209/(300304/(301401/(302500/(303601/
(304704/(305809/(306916/(308025/(309136/(310249/(311364/(312481/(313600/(314721/(315844/(316969/
(318096/(319225/(320356/(321489/(322624/(323761/(324900/(326041/(327184/(328329/(329476/(330625/
(331776/(332929/(334084/(335241/(336400/(337561/(338724/(339889/(341056/(342225/(343396/(344569/
(345744/(346921/(348100/(349281/(350464/(351649/(352836/(354025/(355216/(356409/(357604/(358801/
(360000/(361201/(362404/(363609/(364816/(366025/(367236/(368449/(369664/(370881/(372100/(373321/
(374544/(375769/(376996/(378225/(379456/(380689/(381924/(383161/(384400/(385641/(386884/(388129/
(389376/(390625/(391876/(393129/(394384/(395641/(396900/(398161/(399424/(400689/(401956/(403225/
(404496/(405769/(407044/(408321/(409600/(410881/(412164/(413449/(414736/(416025/(417316/(418609/
(419904/(421201/(422500/(423801/(425104/(426409/(427716/(429025/(430336/(431649/(432964/(434281/
(435600/(436921/(438244/(439569/(440896/(442225/(443556/
(444889/835+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)​+2672)+668)+2672
)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+6​68)+2672)+668)+2672)+668
)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2​672)+668)+2672)+668)+267
2)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+​668)+2672)+668)+2672)+66
8)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+​2672)+668)+2672)+668)+26
72)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)​+668)+2672)+668)+2672)+6
68)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)​+2672)+668)+2672)+668)+2
672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672​)+668)+2672)+668)+2672)+
668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668​)+2672)+668)+2672)+668)+
2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+267​2)+668)+2672)+668)+2672)
+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+66​8)+2672)+668)+2672)+668)
+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+26​72)+668)+2672)+668)+2672
)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+6​68)+2672)+668)+2672)+668
)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2​672)+668)+2672)+668)+267
2)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+​668)+2672)+668)+2672)+66
8)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+​2672)+668)+2672)+668)+26
72)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)​+668)+2672)+668)+2672)+6
68)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)​+2672)+668)+2672)+668)+2
672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672​)+668)+2672)+668)+2672)+
668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668​)+2672)+668)+2672)+668)+
2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+267​2)+668)+2672)+668)+2672)
+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+66​8)+2672)+668)+2672)+668)
+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+26​72)+668)+2672)+668)+2672
)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+6​68)+2672)+668)+2672)+668
)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2​672)+668)+2672)+668)+267
2)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+​668)+2672)+668)+2672)+66
8)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+​2672)+668)+2672)+668)+26
72)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)​+668)+2672)+668)+2672)+6
68)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)​+2672)+668)+2672)+668)+2
672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672​)+668)+2672)+668)+2672)+
668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668​)+2672)+668)+2672)+668)+
2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+267​2)+668)+2672)+668)+2672)
+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+66​8)+2672)+668)+2672)+668)
+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+26​72)+668)+2672)+668)+2672
)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+668)+2672)+6​68)+2672)+668)+2672)+668
)+2672)+668)+2672)+668)+2672)+668)+2672))*4'

EXPAND

-->


209825968201627626286623454595828256626709353294427565693247490908723674963328
437446737220228925859354323920568231930674733184239295858288420011592440669384
872031012213258537332909445846000903747665658888180242434346594380978304852173
052403372808186145692544868234515711423406188811364772876865637765490919076434
392771210399189504646019426036882136493770586097156188551653562529227766165281
757151724133041036460071532563393677368044337870650100262757747252991375882380
250620201017085686160848392397092831545551473286217654554917628548452574061153
069975606948031055019988881584904506395644140050798077890023084451291894497326
051903961793213476739142633750444261417276165846644921102222993269581502224908
997341571378437977195178421097830535179345748645267789093772601974921499991845
401937587852170647321448405827097658859639498298473057686922497943250306778251
041916305052893735026021428385807394291381395492996962793986222406725018359953
072078004002223662842307399263336028649795843721429601453448755585367626557498
760685012488702797279250113419102483986306543454752334763491632369843706569991
047137574804721539328768935056574597269592669763395749629148319846053893268546
966267252153515482975568525561183612983615350523775562374113789568264561165733
980089432434864965104974017069056/
667896800566637703395851118513503626141262321582626289553702325261623907354572
438254043502004959002868854942975857307918561401611054958654114240264009515552
200923163423327214742168578084961721459667457020661333684638926852020031848136
506591215222804140858795351641744040655966457807513342355269848746467904921756
092395537006240650450942873367057350825035314920399047941757674877322573778408
806778939307044804124356789509459360093717241638776925160070406891674464351926
499274636660943669592601023521686228225528199716726528231684826695828870347304
536961944452922511457510213498499968525892553651189414843449897338673657505718
626141672674142307483362631546562482632710881847919379078547871416496195345019
159779507094766464828561525894215976515907672692661978234580958904825570303785
429435771125352380855627325739133377789931371074249342603954029828721762250029
731123836036730436729508866487381850274816327061567345996978616965553971425079
225277549473507787885520688837580635372146777114593848029172095673038532289136
146665525606818535828272746221383566300894538111563291705413589743201178836426
902316517664023633837107436210497878988982635420985619424857564520193944361919
714122045000936522301890676982272649077234522014303978345842413157035790893429
82462406016044012038707143579375


1000 SWAP
FXND DUP SIZE R->I ALOG OVER - PICK3 * SWAP IQUOT + ->STR DUP HEAD -51
FC? { "." } { "," } IFTE + SWAP TAIL + 1 ROT 2 + SUB

EVAL

-->

3.1415926535897932384626433832795028841971693993751058209749445923078164062862
089986280348253421170679821480865132823066470938446095505822317253594081284811
174502841027019385211055596446229489549303819644288109756659334461284756482337
867831652712019091456485669234603486104543266482133936072602491412737245870066
063155881748815209209628292540917153643678925903600113305305488204665213841469
519415116094330572703657595919530921861173819326117931051185480744623799627495
673518857527248912279381830119491298336733624406566430860213949463952247371907
021798609437027705392171762931767523846748184676694051320005681271452635608277
857713427577896091736371787214684409012249534301465495853710507922796892589235
420199561121290219608640344181598136297747713099605187072113499999983729780499
510597317328160963185950244594553469083026425223082533446850352619311881710100
031378387528865875332083814206171776691473035982534904287554687311595628638823
537875937519577818577805321712268066130019278766111959092164201989

Edited to (try to) fix some formatting

[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

[Gerson W. Barbosa]

(02-15-2018 03:00 AM)Mike (Stgt) Wrote:  

(02-14-2018 11:57 PM)Gerson W. Barbosa Wrote:  I believe the continued fraction which works as a correction to the very slowly convergent Madhava-Leibniz series (a.k.a. Gregory series), was known to Madhava.

If you complain about the "very slowly convergent Madhava-Leibniz series" then give the Nilakantha's Series a try... looks nice, but only theoretically convergent, but not realy. (Well, I still have to grasp this paper.)

Thanks for the link! (even though it is surely more difficult for me to grasp it than to you).

Although the plain Madhava-Leibnitz series is very slowly convergent to be of any practical use, when properly tweaked it can be used to compute one thousand digits or more on a relatively slow calculator like a 49G, for instance, in a reasonable time (by reasonable time I mean before the batteries drain out :-). I have first become aware of this fact by one of the Valentin Albillo's legendary mini-challenges. I think I could take that one now, but I fear it's a bit too late :-)

Madhava himself gave three correction formulas to be used after a certain number of terms are evaluated. As I've already said, it is quite possible that he could have given as many of these correction terms as he wanted, by simply calculating approximants to the continued fraction. I found it empirically, by trial and error, but this can be done systematically, as show in this paper from Research India Publications.

Gerson.

[Gerson W. Barbosa]

(02-16-2018 02:32 AM)Mike (Stgt) Wrote:  

(02-15-2018 11:57 PM)Gerson W. Barbosa Wrote:  Algorithm:

 d = 1000                                       ; desired number of digits
 n = int((d + 4)/3)                             ; required number of iterations
...
This gives \(3\) digits per iteration.

Sorry, under what system did you run this algorithm? I did edit/copy-edit/paste it to REXX, few minor adjustments (Int() replaced by integer devision by 1, For/Next loop replaced by Do/End - BTW, the loop counter i is not used, so a Do n without i = 1 to is sufficient, and the command to bring the result to the CRT differs).
Result far off from Pi. So I tried with 10 digits: 0.6170723428 what is much closer but still not fully satisfactory.

Ciao.....Mike

Sorry! There’s a typo in the first line after the end of loop. It should read

s2 = 1/(b + s2)

I will fix it.

Gerson.

PS.: I’ve just tested it on the 75, with n = int(d/3) + 1

1: 3.14285714286
3: 3.1415935416
6: 3.1415926543
9: 3.14159265359

[Gerson W. Barbosa]

(02-16-2018 03:39 AM)Mike (Stgt) Wrote:  

(02-16-2018 03:04 AM)Gerson W. Barbosa Wrote:  s2 = 1/(b + s2)

Works. (Was that a test who's following the thread?)

Ciao.....Mike

That’s essentially the same algorithm in my second program in post #38. For some reason I decided to rename S0 and S1 to S1 and S2, respectively, but my doublecheck (or lack of it) has failed. Hopefully the RPL program is ok (just a copy & paste matter).

Gerson

[EdS2]

I just came across this nice approximation, by Ramanujan (of course)
∜(2143/22) = 3.14159265258...
(more here including some different ways to express the same result.)

Edit: this same approximation is mentioned in these older threads:
CALCULATING MANY DIGITS OF PI (2004)
Computing PI -- a simple programming problem (2006)
Yet another 12C mini-challenge (pi) (2009)

[Gerson W. Barbosa]

(02-17-2018 12:27 PM)EdS2 Wrote:  I just came across this nice approximation, by Ramanujan (of course)
∜(2143/22) = 3.14159265258...

That's probably the only Ramanujan approximation that doesn't rely on any of his astounding theories. He just noticed that \(\pi ^{4}\) = 97.0409091034..., very close to 97.0409090909..., which when multiplied by 990 this gives 96435. Thus, \(\pi \approx \sqrt[4]{\frac{96435}{990}}\), or \(\pi \approx \sqrt[4]{\frac{2143}{22}}\).

Likewise, \((e^{\pi })^{4}\) = \(e^{4\pi }\) = 286751.31313665..., which when multiplied by 99 gives 28388380.0005287... However, \(\frac{28388380}{99}\) is irreducible. Still, \(\frac{\ln \left ( \frac{28388380}{99} \right )}{4 }\) = 3.141592653585137 is an approximation to \(\pi\), although not nearly as good as Ramanujan's, as 10 digits and two operations are used to produce only 12 digits.

For \(e^{\pi }\) I would suggest these two approximations:

\(\frac{16\ln 878}{\ln\left ( 16\ln 878 \right )}\) = 23.14069263691337

and

\(\frac{64146}{2772}\) = 23.14069264069264

the latter being a palindromic approximation.

Of course, the natural logarithm of these are also approximations to \(\pi\).

Gerson.

Edited to fix a typo.

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

PS:

Like the ones for \(\pi\) and for \(e ^{\pi }\), the fourth power of \(\ln \pi\), 1.71716522553..., will make for another approximation. So now we have


\(\pi \approx \sqrt[4]{\frac{96435}{990}}\), or \(\pi \approx \sqrt[4]{\frac{2143}{22}}\) = 3.141592652582646 (3589793)

\(e^{\pi }\approx \sqrt[4]{\frac{28388380}{99}}\) = 23.14069263267153 (277926)

\(\ln \pi \approx \sqrt[4]{\frac{170}{99}}\) = 1.14473096774 (2988585)

The latter can be improved to obtain yet another approximation for \(\pi\), albeit a not so good one:

\(\pi \approx e ^{\sqrt[4]{\frac{170-\frac{400}{622401}}{99}}}\) = 3.1415926535897932121 (384)