HP Forums
Happy Pi Day! (again) - Printable Version

+- HP Forums (https://archived.hpcalc.org/museumforum)
+-- Forum: HP Museum Forums (https://archived.hpcalc.org/museumforum/forum-1.html)
+--- Forum: Old HP Forum Archives (https://archived.hpcalc.org/museumforum/forum-2.html)
+--- Thread: Happy Pi Day! (again) (/thread-188443.html)



Happy Pi Day! (again) - Gerson W. Barbosa - 07-22-2011

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) - Egan Ford - 07-22-2011

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) - Gerson W. Barbosa - 07-22-2011

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) - Oliver Unter Ecker - 07-22-2011

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) - Gerson W. Barbosa - 07-22-2011

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


Re: Happy Pi Day! (again) - Egan Ford - 07-23-2011

Quote:
Notice that the post time is 3:14 pm ...

I missed that. Very clever.


Re: Happy Pi Day! (again) - Gerson W. Barbosa - 07-23-2011

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) - Gerson W. Barbosa - 07-24-2011

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) - Gerson W. Barbosa - 07-29-2011

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.