%I #97 Dec 21 2018 00:07:35
%S 0,1,3,5,6,8,10,11,13,15,16,18,20,21,23,25,26,28,29,31,33,34,36,38,39,
%T 41,43,44,46,48,49,51,53,54,56,58,59,61,63,64,66,68,69,71,73,74,76,78,
%U 79,81,83,84,86,88,89,91,93,94,96,98,99,101,103,104,106,108,109,111,113,114
%N Consider the method used by Archimedes to determine the value of Pi (A000796). This sequence denotes the number of iterations of his algorithm which would result in a difference of less than 1/10^n from that of Pi.
%C It takes on average 5/3 iterations to yield another digit in the decimal expansion of Pi.
%C The side of a 96-gon inscribed in a unit circle is equal to sqrt(2-sqrt(2+sqrt(2+sqrt(2+sqrt(3))))). This is the size of one of the two polygons that Archimedes used to derive that 3 + 10/70 < Pi < 3 + 10/71.
%C In the Mathematica program, we started with an inscribed triangle and a circumscribed triangle of a unit circle and used decimal precision to just over a 1000 places.
%C The perimeter of the circumscribed 3*2^n-polygon exceeds Pi by more than the deficit of the perimeter of the inscribed 3*2^n-polygon. If we were to give twice the weight of the inscribed 3*2^n-polygon to that of the circumscribed 3*2^n-polygon, then the convergence would be twice as fast!
%C From _A.H.M. Smeets_, Jul 12 2018: (Start)
%C Archimedes's scheme: set upp(0) = 2*sqrt(3), low(0) = 3 (hexagons); upp(n+1) = 2*upp(n)*low(n)/(upp(n)+low(n)) (harmonic mean, i.e., 1/upp(n+1) = (1/upp(n) + 1/low(n))/2), low(n+1) = sqrt(upp(n+1)*low(n)) (geometric mean, i.e., log(low(n+1)) = (log(upp(n+1)) + log(low(n)))/2), for n >= 0. Invariant: low(n) < Pi < upp(n); variant function: upp(n)-low(n) tends to zero for n -> inf. The error of low(n) and upp(n) decreases by a factor of approximately 4 each iteration, i.e., approximately 2 bits are gained by each iteration.
%C From Archimedes's scheme, set r(n) = (2*low(n) + upp(n))/3, r(n) > Pi and the error decreases by a factor of approximately 16 for each iteration, i.e., approximately 4 bits are gained by each iteration. This is often called "Snellius acceleration".
%C For similar schemes see also A014549 (in this case with quadratically convergence), A093954, A129187, A129200, A188615, A195621, A202541.
%C Note that replacing "5/3" by "log(10)/log(4)" would be better in the first comment. (End)
%D Petr Beckmann, A History of Pi, 5th Ed. Boulder, Colorado: The Golem Press (1982).
%D Jonathan Borwein and David Bailey, Mathematics by Experiment, Second Edition, A. K. Peters Ltd., Wellesley, Massachusetts 2008.
%D Jonathan Borwein & Keith Devlin, The Computer As Crucible, An Introduction To Experimental Mathematics, A. K. Peters, Ltd., Wellesley, MA, Chapter 7, 'Calculating [Pi]' pp. 71-79, 2009.
%D Eli Maor, The Pythagorean Theorem, Princeton Science Library, Table 4.1, page 55.
%D Daniel Zwillinger, Editor-in-Chief, CRC Standard Mathematical Tables and Formulae, 31st Edition, Chapman & Hall/CRC, Boca Raton, London, New York & Washington, D.C., 2003, ยง4.5 Polygons, page 324.
%H Mike Bertrand, Ex Libris, <a href="http://nonagon.org/ExLibris/archimedes-pi"> Archimedes and Pi</a>
%H Frits Beukers and Weia Reinboud, <a href="https://dspace.library.uu.nl/handle/1874/26396">Snellius versneld</a>, (text in English), preprint.
%H Frits Beukers and Weia Reinboud, <a href="http://www.nieuwarchief.nl/serie5/pdf/naw5-2002-03-1-060.pdf">Snellius versneld</a>, (text in English), NAW 5/3 no. 1, pp. 60-63 (2002).
%H Lee Fook Loong Eugene, <a href="http://www.math.nus.edu.sg/~matlhh/UROPS/reports/LFLreport.pdf">The Computation of [Pi] And Its History</a>
%H Kyutae Paul Han, <a href="https://math.dartmouth.edu/~m56s13/Han_proj.pdf"> Pi and Archimedes Polygon Method</a>
%H Eric Weisstein's World of Mathematics, <a href="http://mathworld.wolfram.com/ArchimedesRecurrenceFormula.html">Archimedes' Recurrence Formula</a>
%H Eric Weisstein's World of Mathematics, <a href="http://mathworld.wolfram.com/RegularPolygon.html">Regular Polygon</a>
%H Michael Woltermann Ph.D., Washington & Jefferson College, <a href="http://www2.washjeff.edu/users/mwoltermann/Dorrie/38.pdf">38. Archimedes' Determination of [Pi].</a>
%F Conjecture: There exists a c such that a(n) = floor(n*log(10)/log(4)+c); where c is in the range [0.08554,0.10264]. Critical values to narrow the range are believed to be at a(74), a(133), a(192), a(251), a(310), a(366), a(425), a(484). - _A.H.M. Smeets_, Jul 23 2018
%e Just averaging the initial two triangles (3.89711) would yield Pi to one place of accuracy, i.e., the single digit '3'. Therefore a(0) = 0.
%e The first iteration yields, as the perimeters of the two hexagons, 4*sqrt(3) and 6. Their average is ~ 3.2320508 which is within 1/10 of the true value of Pi. Therefore a(1) = 1.
%e a(3) = 5 since it takes 5 iterations of Archimedes's algorithm to drive the averaged value of the circumscribed 96-gon and the inscribed 96-gon to yield a value within 0.001 of the correct value of Pi.
%e a(4) = 6 since it takes 6 iterations of Archimedes's algorithm to drive the averaged value of the circumscribed 3*2^6-gon and the inscribed 3*2^6-gon to yield a value within 0.0001 of the correct value of Pi.
%t a[n_] := a[n] = N[2 a[n - 1] b[n - 1]/(a[n - 1] + b[n - 1]), 2^10]; b[n_] := b[n] = N[ Sqrt[ b[n - 1] a[n]], 2^10]; a[-1] = 2Sqrt[27]; b[-1] = a[-1]/2; f[n_] := Block[{k = 0}, While[ 10^n*((a[k] + b[k])/4 -Pi) > 1, k++]; k]; Array[f, 70]
%Y Cf. A000796.
%K nonn,base,easy
%O 0,3
%A _William H. Richardson_ and _Robert G. Wilson v_, Jul 03 2014