The Museum of HP Calculators

HP Forum Archive 20

[ Return to Index | Top of Index ]

Happy Pi Day! (again)
Message #1 Posted by Gerson W. Barbosa on 22 July 2011, 3:14 p.m.

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

      
Re: Happy Pi Day! (again)
Message #2 Posted by Egan Ford on 22 July 2011, 3:26 p.m.,
in response to message #1 by Gerson W. Barbosa

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.

            
Re: Happy Pi Day! (again)
Message #3 Posted by Gerson W. Barbosa on 22 July 2011, 3:43 p.m.,
in response to message #2 by Egan Ford

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:

http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/archv019.cgi?read=146977

Too bad the picture is missing. It was a nice heart!

Regards,

Gerson.

                  
Re: Happy Pi Day! (again)
Message #4 Posted by Egan Ford on 23 July 2011, 2:42 p.m.,
in response to message #3 by Gerson W. Barbosa

Quote:
Notice that the post time is 3:14 pm ...
I missed that. Very clever.
                        
Re: Happy Pi Day! (again)
Message #5 Posted by Gerson W. Barbosa on 23 July 2011, 3:06 p.m.,
in response to message #4 by Egan Ford

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.

      
Re: Happy Pi Day! (again)
Message #6 Posted by Oliver Unter Ecker on 22 July 2011, 6:55 p.m.,
in response to message #1 by Gerson W. Barbosa

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.
            
Re: Happy Pi Day! (again)
Message #7 Posted by Gerson W. Barbosa on 22 July 2011, 7:24 p.m.,
in response to message #6 by Oliver Unter Ecker

...using Vieta's formula, if I were to guess.

      
Re: Happy Pi Day! (again)
Message #8 Posted by Gerson W. Barbosa on 24 July 2011, 11:03 p.m.,
in response to message #1 by Gerson W. Barbosa

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.
            
Re: Happy Pi Day! (again)
Message #9 Posted by Gerson W. Barbosa on 29 July 2011, 8:00 p.m.,
in response to message #8 by Gerson W. Barbosa

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.


[ Return to Index | Top of Index ]

Go back to the main exhibit hall