Calculating pi

Here’s a lovely iteration which converges quadratically to pi. Start with:

x_0=sqrt{2}

pi_0=2+sqrt{2}

y_1=2^{1/4}.

Then:

displaystyle{x_{n+1}=frac{1}{2}left(sqrt{x_n}+frac{1}{sqrt{x_n}}right)} for nge 0

displaystyle{y_{n+1}=frac{y_nsqrt{x_n}+1/sqrt{x_n}}{y_n+1}} for nge 1

displaystyle{pi_n=pi_{n-1}frac{x_n+1}{y_n+1}} for nge 1

This remarkable iteration comes from Borwein and Borwein’s “Pi and the AGM“; this is the first formula for pi of the many in the book. It is based on the arithmetic-geometric mean, which is defined as follows:

Given a,b define the sequences a_n,b_n by a_0=a,b_0=b and for all nge 1

displaystyle{a_n=frac{a_{n-1}+b_{n-1}}{2}}

displaystyle{b_n=sqrt{a_{n-1}b_{n-1}}}.

It can be shown that these sequences both converge to the same value, called the arithmetic-geometric mean of a and b, and denoted

M(a,b).

The derivation of the iteration above starts with the complete elliptic integral of the first kind:

displaystyle{K(k)=int^{pi/2}_0frac{dtheta}{sqrt{1-k^2sin^2theta}}}

The importance of these integrals for the AGM includes the useful result

M(1,k)=pi K(sqrt{1-k^2})/2.

It can also be shown that:

displaystyle{Kleft(frac{1}{sqrt{2}}right)dot{K}left(frac{1}{sqrt{2}}right)=frac{pi}{sqrt{2}}}

and that

displaystyle{dot{M}(1,k)=frac{pi}{2}frac{k}{sqrt{1-k^2}}frac{dot{K}(sqrt{1-k^2})}{K^2(sqrt{1-k^2})}}

(In these results, the upper dot indicates the derivative. It is conventional to write k' for sqrt{1-k^2} but to keep the symbols down I won’t use that here.) Finally, and this is the basis for the iteration:

displaystyle{pi=2sqrt{2}frac{M^3(1,1/sqrt{2})}{dot{M}(1,1/sqrt{2})}}.

Full details are given in the book above.

Now, let’s see this in operation. Using Maxima, with floating point precision set at 200 digits; we start with:

(%1) fpprec:200;
(%i2) [x,y,pi]:bfloat([sqrt(2),2^(1/4),2+sqrt(2)])pi; (%o2) 3.4142135623730950488016887242096980785696718753769 48073176679737990732478462107038850387534327641572735013846 23091229702492483605585073721264412149709993583141322266592 75055927557999505011527820605715b0 </code>  Now lets apply the iteration a few times, each time checking the number of digits againstlatex pi:  <code> (%i3) for i:1 thru 6 do (x:(sqrt(x)+1/sqrt(x))/2, y1:(y*sqrt(x)+1/sqrt(x))/(y+1), pi:pi*(x+1)/(y+1),print("pi = ",pi),print("diff = ", floor(abs(log(pi-bfloat(%pi))/log(10.0)))),print(),y:y1); </code> <code>pi =  3.142606753941622600790719823618301891971356246277167 25391106706733002432829841433917490844363218215939832719424 63598401856831575329982467730414626149809288536376114352367 443809738134706288502949125099b0 diff =  2 </code> <code>pi =  3.141592660966044230497752235120339690679284256864528 90583358376281661542951772210269832001264427102653464852593 19110073690752404626908089759117363527617678835997985181167 4320514766167616002203773916b0 diff =  8 </code> <code>pi =  3.141592653589793238645773991757141794034789623867451 84194317618340870893816338362721980735705521698725721871080 51851694141568019749840320178331947918544614132145382957429 175251304327144303796308774166b0 diff =  18 </code> <code>pi =  3.141592653589793238462643383279502884197224120466562 72039327207763960604807177493746574457598832522036441773314 43692556945240355112240152544528702433558739837465363838916 511954888511539934938771670765b0 diff =  40 </code> <code>pi =  3.141592653589793238462643383279502884197169399375105 82097494459230781640628620899863044094747256210005455173392 24142307579850000984756761462246969993739564082271697581825 398673085036296812176348099349b0 diff =  83 </code> <code>pi =  3.141592653589793238462643383279502884197169399375105 82097494459230781640628620899862803482534211706798214808651 32823066470938446095505822317253594081284811174502841027019 408296862745790132157762511163b0 diff =  170 </code> In only six iterations we have 170 correct digits.  Since the number of digits approximately doubles at each step, if we could perform arbitrarily precise computations, we could produce a million digits oflatex piin only 13 more iterations.  Borwein and Borwein give many other formulas and iterations forlatex pi$, including one which exhibits quintic convergence; that is, the number of correct digits multiplies by five at each step. To find that you'll have to read their book...

2 thoughts on “Calculating pi

  1. Possible typo?

    > the number of correct digits multiples by five at each step.
    > the number of correct digits multiplies by five at each step.

  2. Wow, I just added Brent-Salamin to my blog … then took a look at yours. Your algorithm looks very similar to mine, though seems simpler?

Leave a Reply

Your email address will not be published. Required fields are marked *