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

.

Then:

for

for

for

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

Given define the sequences by and for all

.

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

.

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

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

.

It can also be shown that:

and that

(In these results, the upper dot indicates the derivative. It is conventional to write for but to keep the symbols down I won’t use that here.) Finally, and this is the basis for the iteration:

.

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:

(%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 :

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

25391106706733002432829841433917490844363218215939832719424

63598401856831575329982467730414626149809288536376114352367

443809738134706288502949125099b0

diff = 2

90583358376281661542951772210269832001264427102653464852593

19110073690752404626908089759117363527617678835997985181167

4320514766167616002203773916b0

diff = 8

84194317618340870893816338362721980735705521698725721871080

51851694141568019749840320178331947918544614132145382957429

175251304327144303796308774166b0

diff = 18

72039327207763960604807177493746574457598832522036441773314

43692556945240355112240152544528702433558739837465363838916

511954888511539934938771670765b0

diff = 40

82097494459230781640628620899863044094747256210005455173392

24142307579850000984756761462246969993739564082271697581825

398673085036296812176348099349b0

diff = 83

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

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

Possible typo?

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

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

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?