Happy Pi Day! (again) « Next Oldest | Next Newest »

 ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 07-22-2011, 03:14 PM Pi Day is a harmless nerdity, so no problem having two of them in the same year, 130 days apart from each other :-) For people using the European date system (dd/mm/yyyy), 22/07 matches exactly Archimedes' upper bound for pi. Without algebra, trigonometry and a positional number system, Archimedes was able to place pi between the bounds 6336/20171/4 and 14688/46731/2. He widened very slightly the interval and established the more simple approximations 310/71 and 31/7, the latter being the de facto value for pi in Europe during more than fifteen centuries, so much that some believed that was an exact value, not an approximation. Archimedes used Euclidean geometry. Essentially, he established those bounds by calculating the perimeters of two 96-gon polygon, one inscribed and the other circumscribed to a one-unit diameter circumference, having started with hexagons and doubling the number of the polygons five times. There are evidences went even further, but because accumulated copying errors along the centuries the bounds in a surviving palimpsest are way off: 211875/67441 and 197888/62351. 211872/67441 and 195888/62353 are plausible reconstructed values and implied Archimedes might have gone up the 1536-gon. The Archimedes' Method was the first analytical method for computing pi. It is far from being efficient, but it is simple and beautiful nonetheless. It offers only linear convergence, giving 0.6 decimal digits per iteration and requires the extraction of square roots. However, there are hidden relationships between the bounds that allows for two-, three- and four-fold speed-ups, perhaps even more. In another post, I have presented some programs and a couple of formulas that give 1.8 digits per iteration. I've found a relationship between them and had WolframAlpha solve the equation (just submit "solve (p-(27*(a^2*b)^(1/3)-4*(2*a+b))/15 )/((9*a*c/(4*a-c)-7/3*(a-4*c))/10-p)==128/19 for p" to W|A). The result can be simplified to: ``` pi ~ (513*(a^2*b)^(1/3) - 24*a*(288*a/(c - 4*a) + 97) - 76*b + 1792*c)/2205 where a = lower bound b = upper bound c = next lower bound ``` This gives 2.4 digits per iteration. For checking purpose only, we can use trigonometry and evaluate the following: ```n:=14 a:=3*2^n*sin(pi/(3*2^n)) b:=3*2^n*tan(pi/(3*2^n)); c:=3*2^(n+1)*sin(pi/(3*2^(n+1))) p:=(513*(a^2*b)^(1/3)-24*a*(288*a/(c-4*a)+97)-76*b+1792*c)/2205``` On W|A the input line "compute {n=14, p=(513*(a^2*b)^(1/3)-24*a*(288*a/(c-4*a)+97)-76*b+1792*c)/2205, a=3*2^n*sin(pi/(3*2^n)), b=3*2^n*tan(pi/(3*2^n)), c=3*2^(n+1)*sin(pi/(3*2^(n+1)))} to 42 places" will return: ```a ~ 3.141592651450767651704253640492219020448 b ~ 3.141592657867844419844008529336759834985 c ~ 3.141592653055036841691123180415474202257 n ~ 14.000000000000000000000000000000000000000 p ~ 3.14159265358979323846264338327950288406 (our underlining emphasizes the matching digits) ``` Notice that the perimeters of the inscribed and circumscribed 49152-gons (3*2^14) and the perimeter of the inscribed 98304-gon give only 9 digits of accuracy and yet the relatively simple expression involving them, evaluated only in the last iteration, gives out pi to 36 significant digits. As a comparison, in early 17th century Ludolph Van Ceulen used a polygon with 2^62 sides to achieve the same accuracy. If he had investigated how the bounds related to each others, for only the first few iterations, he might have discovered at least the simple weighted mean, which would have allowed him to find the first 50 digits. The lack of mechanical calculators, however, did not encourage numerical investigations. The Excel sheet excerpt below may require some editing for clarity, but it gives an idea how the limits can be found numerically: ``` A B C D E F G H I J 1 (lower bound) (upper bound) 2 n 3*2^n (3*2^n*sin(pi/(3*2^n)) (3*2^n*tan(pi/(3*2^n)) (pi - an)/(pi - an+1) (pi/an - 1)/(pi/an+1 - 1) (4*an+1 - an)/3 (3*an*an+1)/(4*an - an+1) (H-pi)/(pi-G) (7*G + 3*H)/10 3 (a) (b) 4 0 3 2,59807621135332 5,19615242270663 3,83859210528741 4,43242437059372 3,13397459621556 3,16311169400545 2,82474118512721 3,14271572555253 5 1 6 3,00000000000000 3,46410161513775 3,95907081843196 4,09873171487929 3,14110472164033 3,14278367587695 2,44095982743118 3,14160840791132 6 2 12 3,10582854123025 3,21539030917347 3,98973131852282 4,02415855279584 3,14156197063157 3,14166504789936 2,35943056841153 3,14159289381191 7 3 24 3,13262861328124 3,15965994209750 3,99743055062314 4,00600772065997 3,14159073296874 3,14159714747608 2,33980893047569 3,14159265732094 8 4 48 3,13935020304687 3,14608621513143 3,99935749514110 4,00149994832716 3,14159253350506 3,14159293398155 2,33494921116754 3,14159265364801 9 5 96 3,14103195089051 3,14271459964537 3,99983936486646 4,00037486339965 3,14159264608378 3,14159267110686 2,33373691378455 3,14159265359070 10 6 192 3,14145247228546 3,14187304997982 3,99995984066625 4,00009370815163 3,14159265312066 3,14159265468449 2,33343335829224 3,14159265358981 11 7 384 3,14155760791186 3,14166274705685 3,99998995986230 4,00002342619481 3,14159265356047 3,14159265365821 2,33330304269466 3,14159265358979 12 8 768 3,14158389214832 3,14161017660469 -> 4 -> 4 -> 7/3 Solve Solve Solve (pi - an)/(pi - an+1) = 4 (pi/an - 1)/(pi/an+1 - 1) = 4 (H-pi)/(pi - G) = 7/3 for pi to get F for pi to get G for pi to get J ``` If only Van Ceulen had a spreadsheet application... Definitely calculating devices were very limited back then :-) Edited to fix some typos Edited: 23 July 2011, 3:18 p.m. after one or more responses were posted ▼ Egan Ford Posting Freak Posts: 1,619 Threads: 147 Joined: May 2006 07-22-2011, 03:26 PM Happy Pi Day to you too! Excellent article. Thanks. BTW, for those that may have missed my "official" Pi Day post here it is again: Pi Day Rematch: Apple II vs. HP-41C and a follow-up: Post Pi Day Post. ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 07-22-2011, 03:43 PM Thank you, Egan! Notice that the post time is 3:14 pm and the edit time 3:18, the digits of 1/pi since today is 7/22 also. Well, the latter was not intentional but the explanation fits nicely :-) I certainly didn't miss your really excellent Pi Day post. I have yet to run Robert Bishop's Apple PI on the emulator. By the way, my first computer was a somewhat limited Apple II clone. When I was studying Calculus I wrote a graph plotting program for it: Too bad the picture is missing. It was a nice heart! Regards, Gerson. ▼ Egan Ford Posting Freak Posts: 1,619 Threads: 147 Joined: May 2006 07-23-2011, 02:42 PM Quote: Notice that the post time is 3:14 pm ... I missed that. Very clever. ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 07-23-2011, 03:06 PM Or just very nerdish :-) Thanks for the follow-up to your Pi Day article. I knew about Shank's mistake but I had never seen it in print. Oliver Unter Ecker Member Posts: 239 Threads: 9 Joined: Dec 2010 07-22-2011, 06:55 PM For your enjoyment: This ``` 6666,-2%{2+.2/@*\/10.3??2*+}* ``` generates the first 1000 digits of Pi in GolfScript, a stack-based language. Taken from GolfScript Examples. ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 07-22-2011, 07:24 PM ...using Vieta's formula, if I were to guess. Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 07-24-2011, 11:03 PM This approach can be carried out one step further when we eventually obtain a formula that yields 3 digits per iteration: ``` pi ~ (640*c*(28*c-139*a)/(c-4*a)+4*b*(3737+14256*b/(d-4*b))+1056*d-2568*a-6453*(a^2*b)^(1/3))/11655 where a: lower bound b: upper bound c: next lower bound d: next upper bound ``` This is a five-fold improvement over the original method. At least two iterations have to be executed since the approximation formula requires two successive sets of bounds. Here are the results obtained from the first four iterations, starting with hexagons and dodecagons: ``` 6-12 3.141592925577355371 12-24 3.141592653812338501 24-48 3.141592653590001584 48-96 3.141592653589793440 pi 3.141592653589793238 ``` The digits per iteration ratio approaches asymptotically 3 as the number of iterations increases. It is a bit greater initially, as we can see. Notice that the 48-gon suffices for the approximation to pi displayed in 12-digit calculators. ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 07-29-2011, 08:00 PM My HP-50g has boldly computed the first 1000 decimal places of pi using this algorithm (up to the 6*2^333-gon, 10.5 times more sides than a googolgon!) in "only" 2h 59m 11s, thanks to the use of the formula above. Otherwise, it would have taken about fifteen hours to go up to the required 6*2^1667-gon. It is a good calculator, I won't ask it to "compute to the last digit of pi" :-) As a comparison, the second program (Brent-Salamin algorithm) does it in only four minutes. Both programs require the LongFloat library for the HP-49G/50g by Gjermund Skailand and Thomas Rast. http://www.hpcalc.org/details.php?id=5363 ```%%HP: T(3)A(R)F(.); DIR QUICKARCH \<< DUP 3 + 'DIGITS' STO 3 / 1 - 3 DUP FSQRT 2 FMULT ROT 1 SWAP START OVER DUP2 FMULT 2 FMULT UNROT FADD FDIV SWAP OVER FMULT FSQRT SWAP NEXT DUP2 OVER DUP2 FMULT 2 FMULT UNROT FADD FDIV SWAP OVER FMULT FSQRT \-> a b d c \<< 28 c FMULT 139 a FMULT FSUB c 4 a FMULT FSUB FDIV 640 FMULT c FMULT b d 4 b FMULT FSUB FDIV 14256 FMULT 3737 FADD 4 b FMULT FMULT FADD 1056 d FMULT FADD 2568 a FMULT FSUB a DUP FMULT b FMULT 3 FXROOT 6453 FMULT FSUB 11655 FDIV ZZ\<-\->F DROP \->STR DUP HEAD "." + SWAP TAIL + DUP SIZE 1 - " " REPL \>> \>> DIGITS 1002 END 1000 QUICKARCH -> 3.1415926535897932384626433832795 028841971693993751058209749445923 078164062862089986280348253421170 679821480865132823066470938446095 505822317253594081284811174502841 027019385211055596446229489549303 819644288109756659334461284756482 337867831652712019091456485669234 603486104543266482133936072602491 412737245870066063155881748815209 209628292540917153643678925903600 113305305488204665213841469519415 116094330572703657595919530921861 173819326117931051185480744623799 627495673518857527248912279381830 119491298336733624406566430860213 949463952247371907021798609437027 705392171762931767523846748184676 694051320005681271452635608277857 713427577896091736371787214684409 012249534301465495853710507922796 892589235420199561121290219608640 344181598136297747713099605187072 113499999983729780499510597317328 160963185950244594553469083026425 223082533446850352619311881710100 031378387528865875332083814206171 776691473035982534904287554687311 595628638823537875937519577818577 805321712268066130019278766111959 092164201989 BSP: %%HP: T(3)A(R)F(.); \<< 1 + 'DIGITS' STO .5 R\<-\->F 1 2 FSQRT 2 FDIV 1 DO DUP2 OVER FADD 2 FDIV ROT ROT FMULT FSQRT OVER 4 ROLL FSUB DUP FMULT 4 ROLL 2 FMULT DUP2 FMULT 6 ROLL SWAP FSUB ROT 5 ROLLD 4 ROLLD SWAP ROT 5 ROLL UNTIL 0 F.EQ END FMULT 2 FMULT ROT FDIV NIP ZZ\<-\->F DROP \->STR DUP HEAD "." + SWAP TAIL + \>> ======================================= Program Quick_Archimedes; var a, b, c, d, p: real; i, m: byte; n: integer; function Sqrt(x: Real): real; var s, t: Real; begin if x<>0 then begin s:=x/2; repeat t:=s; s:=(s+x/s)/2 until s=t; Sqrt:=s end else Sqrt:=0 end; function Cbrt(x: Real): real; var c, t: Real; begin if x<>0 then begin c:=x/2; repeat t:=c; c:=(2*c+x/(c*c))/3 until Abs(c-t)<1e-16; Cbrt:=c end else Cbrt:=0 end; begin ReadLn(m); n:=6; a:=3; b:=2*Sqrt(3); ClrScr; WriteLn(' n-gon',' ':10,'a',' ':21,'b',' ':21,'p',#10); for i:=1 to m-1 do begin WriteLn(n:4,a:22:17,b:22:17); b:=2*a*b/(a+b); a:=Sqrt(a*b); n:=2*n end; WriteLn(n:4,a:22:17,b:22:17); d:=2*a*b/(a+b); c:=Sqrt(a*d); p:=(640*c*(28*c-139*a)/(c-4*a)+4*b*(3737+14256*b/(d-4*b))+1056*d-2568*a-6453*Cbrt(a*a*b))/11655; WriteLn(2*n:4,c:22:17,d:22:17,p:22:17) end. n-gon a b p 6 3.00000000000000000 3.46410161513775460 12 3.10582854123024915 3.21539030917347248 24 3.13262861328123820 3.15965994209750049 48 3.13935020304686722 3.14608621513143498 96 3.14103195089050965 3.14271459964536831 3.14159265358979344 > ``` ======================================= Update: FINV FY\|^X has been replaced with FXROOT in the first program. 3 FINV FY\|^X was taking about 30 minutes; 3 FXROOT takes less than 27 seconds. Edited: 4 Aug 2011, 11:43 a.m.

 Possibly Related Threads... Thread Author Replies Views Last Post An amazing day: Giving a talk at HP about their calculators Geir Isene 9 2,981 12-16-2013, 06:14 PM Last Post: aurelio HP Prime: Quite Happy Helge Gabert 13 3,578 12-10-2013, 08:19 AM Last Post: Eddie W. Shore [OT] Mathematica free for Raspberry PI BruceH 32 5,706 11-23-2013, 05:24 AM Last Post: Nick_S Computing pi with the PC-1300S Kiyoshi Akima 0 692 11-17-2013, 12:24 AM Last Post: Kiyoshi Akima Calculating Pi LHH 9 1,956 09-27-2013, 10:50 PM Last Post: Gerson W. Barbosa HHC 2013 Day 2 Highlights Eddie W. Shore 6 1,684 09-23-2013, 04:03 PM Last Post: Kimberly Thompson HHC 2013: Day 1 Highlights Eddie W. Shore 28 5,221 09-23-2013, 03:22 PM Last Post: Brad Barton Visualization of pi Bruce Bergman 13 2,561 08-17-2013, 05:00 PM Last Post: Howard Owen Happy Mother's Day! Eddie W. Shore 1 764 05-12-2013, 11:35 AM Last Post: Walter B happy fibonacci day 5/8/13 Allen 8 2,064 05-09-2013, 01:48 AM Last Post: Gerson W. Barbosa

Forum Jump: