¼

The constant p ¼

4  Pi and the Arithmetic-Geometric Mean

4.1  The AGM sequence

Studied by Lagrange, Gauss and Legendre the Arithmetic-Geometric Mean (or AGM) of two positive numbers a and b is defined to be the common limit of the two sequences { an,bn} n ³ 0 determined by:

a0 = a
,
       b0 = b
an+1 = (an+bn)
2
,
bn+1 =
Ö
 

anbn
 

We assume a ³ b. From the immediate relation

(
Ö
 

a
 
-
Ö
 

b
 
)2 ³ 0    Þ    (a+b)
2
³
Ö
 

ab
 

we have an ³ bn and then an ³ (an+bn)/2 = an+1 and bn+1 = Ö{anbn} ³ bn for every n. Therefore, by induction:

b
£
b1 £ ... £ bn £ bn+1 £ ... £ an+1 £ an £ ... £ a1 £ a
an+1-bn+1
£
an+1-bn = (an-bn)
2
£ (a0-b0)
2n
an+1-bn+1
=
1
2
((an+bn)-2
Ö
 

anbn
 
) = 1
2
(
Ö
 

an
 
-
Ö
 

bn
 
)2 = (an-bn)2
2(
Ö
 

an
 
+
Ö
 

bn
 
)2

From the first two relations we observe that an and bn converge monotonically to a common limit determined by the initial conditions (a,b). We denote the common limit by:

M(a,b) =
lim
n® ¥ 
an =
lim
n® ¥ 
bn

More, from the third relation, the convergence is of order 2 (or quadratic):

ê
ê
ê
an+1-M(a,b)
(an-M(a,b))2
ê
ê
ê
= O(1)

The Arithmetic-Geometric Mean has obvious properties:

M(a,a)
=
a
M(a,b)
=
M( a+b
2
,   __
Öab
 
)
M(la,lb)
=
lM(a,b)
M(1,x)
=
1+x
2
M(1, 2Öx
1+x
)

Example 1 M(Ö2,1) (this value is related to the lemniscate's constant and the following numerical computations were made by Gauss in 1799)

a0
=
1.414213562373095048802,b0 = 1
a1
=
1.207106781186547524401,b1 = 1.189207115002721066717
a2
=
1.198156948094634295559,b2 = 1.198123521493120122607
a3
=
1.198140234793877209083,b3 = 1.198140234677307205798
a4
=
1.198140234735592207441,b4 = 1.198140234735592207439

The rate of convergence is impressive, | a6-b6| < 10-86!!

For more details on the AGM, and for it's generalization on complex numbers, read the very nice article from Cox ([3]).

4.2  The AGM and the complete elliptic integral of the first kind

The main result is the following theorem, where I(a,b) denotes the first kind complete elliptic integral:

Theorem 1 We suppose that a ³ b are positive numbers then:

I(a,b) = ó
õ
p/2

0 
(a2cos2q+b2sin2q)-1/2dq = p
2M(a,b)

Proof : One proof may be given by two substitutions. The first one, u = btanq is classic and transforms the elliptic integral into

I(a,b) = 1
2
ó
õ
+¥

-¥ 
((a2+u2)(b2+u2))-1/2du.

The second substitution v = (u-ab/u)/2 gives (with some care at u = 0):

I(a,b)
=
1
2
ó
õ
+¥

-¥ 
((a2+u2)(b2+u2))-1/2du = 1
2
ó
õ
+¥

-¥ 
2dv

Ö

((a+b)2+4v2)(ab+v2)
       so:
I(a,b)
=
I( a+b
2
,   __
Öab
 
) = I(a1,b1) = ... = I(an,bn) = ... = I(M(a,b),M(a,b))

and the last integral is trivial so:

I(a,b) = I(M(a,b),M(a,b)) = ó
õ
p/2

0 
dq
M(a,b)
= p
2M(a,b)

This theorem is most important and allows to compute the first kind complete elliptic integral with the Arithmetic-Geometric Mean.

4.3  The Legendre relation

The second kind complete elliptic integral is defined as the length of a quarter of an ellipse with major axis a and minor axis b:

J(a,b) = ó
õ
p/2

0 

Ö
 

a2cos2q+b2sin2q
 
dq

If (b,b¢) are linked by the relation b2+b¢2 = 1, we have the famous Legendre's relation:

Theorem 2

L(b) = I(1,b)J(1,b¢)+I(1,b¢)J(1,b)-I(1,b)I(1,b¢) = p
2

Rather technical, it can be proved, like Legendre did, by derivation of the function L(b) with respect to the variable b, the result is the null function. The constant value of this function is given with b = 0 and b¢ = 1.

Theorem 3 Like for the function I(a,b) we may prove :

J(a,b) = 2J æ
ç
è
a+b
2
,   __
Öab
 
ö
÷
ø
-abI æ
ç
è
a+b
2
,   __
Öab
 
ö
÷
ø

Using the AGM iteration's notation we may write the previous relation as:

J(an,bn) = 2J(an+1,bn+1)-anbnI(an+1,bn+1) = 2J(an+1,bn+1)-anbnI(a0,b0)

with some basic manipulations we find:

J(a0,b0) = æ
ç
è
a02- 1
2
¥
å
k = 0 
2k(ak2-bk2) ö
÷
ø
I(a0,b0).

4.4  Expression for p

We apply the Legendre's relation with b = b¢ = 1/Ö2 so:

2I(1,1/Ö2)J(1,1/Ö2)-I(1,1/Ö2)2
=
p
2
æ
ç
è
1- 1
2
¥
å
k = 0 
2k(ak2-bk2) ö
÷
ø
I(1,1/Ö2)
=
J(1,1/Ö2)

one can now eliminate J(1,1/Ö2) in the Legendre's relation:

æ
è
1- ¥
å
k = 0 
2k(ak2-bk2) ö
ø
I(1,1/Ö2)2 = p
2

and from I(a,b) = p/(2M(a,b)) we obtain:

æ
è
1- ¥
å
k = 0 
2k(ak2-bk2) ö
ø
p2
4M(1,1/ Ö2)2
= p
2

we have just proved the very important result:

Theorem 4 With a0 = 1, b0 = 1/Ö2 and { ak,bk} k ³ 0 defined by the arithmetic-geometric mean:

p = 2M(1,1/Ö2)2
1- ¥
å
k = 0 
2k(ak2-bk2)

This formula for p was discovered in 1976 by Eugene Salamin and also at the same time by Richard P. Brent ([1], [2]). The convergence is of order 2 (the number of decimal places is doubled at each iteration) and if one is able to compute M(1,1/Ö2) with efficient algorithm, this leads to a very nice way to compute many digits of p.

4.5  Description of an algorithm using Brent-Salamin formula

With the previous notation the algorithm starts with the initial values:

a = 1,b = 1/Ö2,t = 1/2,k = 1

with (a,b,d,s,t) being multi-precision numbers and k an integer. Then the following calculations are made until the required precision is reached (a-b < e):

s
=
1
2
(a+b)
d
=
a-s
d
=
d2
a
=
s
s
=
s2
t
=
t-2kd
b
=
  ___
Ös-d
 
k
=
k+1

The value for p is then:

p » (a+b)2
2t

In this algorithm, multi-precision squaring, square roots and one division are required. During the last years, many record of computation were achieved with such formulae. Jointly other formulae are used to check the result (for instance the Borwein quartically algorithm is used to compare the result with Brent-Salamin algorithm).

4.6  The Borwein quartic iteration

In [4], the following algorithm is given

yk+1
=
1-(1-yk4)1/4
1+(1-yk4)1/4
ak+1
=
(1+yk+1)4ak-22k+3yk+1(1+yk+1+yk+12),

starting with

y0 = Ö2-1,a0 = 6-4Ö2,

then ak converge quartically to 1/p. This algorithm provides a very fast and efficient way to compute p and is often used with the Brent-Salamin algorithm to establish new record of computation. Of course it's necessary to compute square roots with high precision and this can be done efficiently thanks to the Fast Fourier Transform.

Here are the first iterates :

1
a1
=
3.1415926(462135422821493444319826957743144372...),
1
a2
=
3.1415926535897932384626433832795028841971(146...),

and 1/a3 is about 170 digits, 1/a4 about 694 digits, ... The number of digits increases by a factor of four at each iteration (see [5]).

In the same book we also find quadratic, cubic, quintic, septic, ... algorithms.

4.7  Some records of computation

1988 Kanada gave 201,326,000 decimal places with the Brent-Salamin formula. Only twenty eight iterations of the Arithmetic-Geometric Mean were needed to achieve such an accuracy. The Borwein iteration was used for the verification.

1996 Kanada gave 17,100,000,000 digits using the same algorithm. The computation was performed on a Hitachi parallel computer.

1999 Kanada gave 68 billions of digits with the same algorithm, then 206 billions digits (Sept. 1999).

References

[1]
R.P. Brent, Fast multiple-Precision evaluation of elementary functions, J. Assoc. Comput. Mach., (1976), vol. 23, p. 242-251

[2]
E. Salamin, Computation of p Using Arithmetic-Geometric Mean, Mathematics of Computation, (1976), vol. 30, p. 565-570

[3]
A. Cox, The Arithmetic-Geometric Mean of Gauss, L'enseignement Mathématique, (1984), vol. 30, p. 275-330

[4]
J.M. Borwein and P.B. Borwein, Pi and the AGM - A study in Analytic Number Theory and Computational Complexity, A Wiley-Interscience Publication, New York, (1987)

[5]
J.M. Borwein and P.B. Borwein, Ramanujan and Pi, Scientific American, (1988), p. 112-117


Back to The constant pi


File translated from TEX by TTH, version 2.32.
On 17 Aug 2000, 21:44.