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/2017^{1}/_{4} and 14688/4673^{1}/_{2}. He widened very slightly the interval and established the more simple approximations 3^{10}/_{71} and 3^{1}/_{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 96gon polygon, one inscribed and the other circumscribed to a oneunit 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 1536gon.
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 fourfold speedups, 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*ac)7/3*(a4*c))/10p)==128/19 for p" to WA). 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/(c4*a)+97)76*b+1792*c)/2205
On WA the input line "compute {n=14, p=(513*(a^2*b)^(1/3)24*a*(288*a/(c4*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 49152gons (3*2^14) and the perimeter of the inscribed 98304gon 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) (Hpi)/(piG) (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 (Hpi)/(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
