Calculating pi

Posted on May 18, 2008 by amca01

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 n\ge 0

\displaystyle{y_{n+1}=\frac{y_n\sqrt{x_n}+1/\sqrt{x_n}}{y_n+1}} for n\ge 1

\displaystyle{\pi_n=\pi_{n-1}\frac{x_n+1}{y_n+1}} for n\ge 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 n\ge 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}_0\frac{d\theta}{\sqrt{1-k^2\sin^2\theta}}}

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{K\left(\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=2\sqrt{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

Now lets apply the iteration a few times, each time checking the number of digits against \pi:


(%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);

pi = 3.142606753941622600790719823618301891971356246277167\
25391106706733002432829841433917490844363218215939832719424\
63598401856831575329982467730414626149809288536376114352367\
443809738134706288502949125099b0
diff = 2

pi = 3.141592660966044230497752235120339690679284256864528\
90583358376281661542951772210269832001264427102653464852593\
19110073690752404626908089759117363527617678835997985181167\
4320514766167616002203773916b0
diff = 8

pi = 3.141592653589793238645773991757141794034789623867451\
84194317618340870893816338362721980735705521698725721871080\
51851694141568019749840320178331947918544614132145382957429\
175251304327144303796308774166b0
diff = 18

pi = 3.141592653589793238462643383279502884197224120466562\
72039327207763960604807177493746574457598832522036441773314\
43692556945240355112240152544528702433558739837465363838916\
511954888511539934938771670765b0
diff = 40

pi = 3.141592653589793238462643383279502884197169399375105\
82097494459230781640628620899863044094747256210005455173392\
24142307579850000984756761462246969993739564082271697581825\
398673085036296812176348099349b0
diff = 83

pi = 3.141592653589793238462643383279502884197169399375105\
82097494459230781640628620899862803482534211706798214808651\
32823066470938446095505822317253594081284811174502841027019\
408296862745790132157762511163b0
diff = 170

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 of \pi in only 13 more iterations.

Borwein and Borwein give many other formulas and iterations for \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…

Filed under: Computation, Maxima

« »

2 Responses

  1. mvngu, on May 18th, 2008 at 9:33 pm Said:

    Possible typo?

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

    Reply
  2. Jake, on May 20th, 2008 at 10:47 pm Said:

    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?

    Reply

Leave a Reply