## 1  Iterative algorithms

Modern iterative algorithms are among the most efficient ways to compute billion of digits of p. The first of those algorithms were discovered during the years 1970. We give without proof a selection and a brief description of some important iterations.

### 1.1  Early iterations

Before the discovery of extremely fast converging modern iterations, a few algorithms were already known but with slow rates of convergence.

#### 1.1.1  Viete

Based on the well known formula discovered by François Viète (1540-1603)

2
p
=   æ
ú
Ö

 12

æ
ú
ú
Ö

1
2
+ 1
2
æ
ú
Ö

 12

æ
ú
ú
Ö

1
2
+ 1
2
æ
ú
Ö

 12 + 12 Ö [1/2]

...

an algorithm may be written as

 xk+1
 =
 xkyk
 yk+1
 =
æ
ú
Ö

 2ykyk+1

and

 x0
 =
 1
 y0
 =
 Ö2

then

 lim k® ¥ xk = p2 .

It converges with a geometrical rate and for example it produces

 2x1
 =
 2.8284...
 2x2
 =
 3.(061...)
 2x3
 =
 3.1(214...)
 2x4
 =
 3.1(365...)
 2x5
 =
 3.14(033...)

#### 1.1.2  Archimedes-Pfaff

Archimedes' method is equivalent to the Pfaff-algorithm published in 1800 (Pfaff was Gauss's teacher).

 xk+1
 =
 2xkykxk+yk
 yk+1
 =
 Ö xk+1yk

starting with

 x0 = 2Ö3,y0 = 3,

then

 lim k® ¥ xk = lim k® ¥ yk = p

again, the rate of convergence is geometric with an error O(4-k) (that is you need roughly about 5 iterations to obtain 3 new digits). It can be accelerated by

 lim k® ¥ zk = lim k® ¥ æç è 13 xk+ 23 yk ö÷ ø = p

which also converges geometrically, but with a faster rate of convergence (O(16-k)). The first iterates of the two sequences are

 y1
 =
 3.1(05828541230249...),z1 = 3.14(2349130544656...)
 y2
 =
 3.1(32628613281238...),z2 = 3.141(639056219992...)
 y3
 =
 3.1(39350203046867...),z3 = 3.14159(5540408389...)
 y4
 =
 3.141(031950890509...),z4 = 3.141592(833808795...)
 y5
 =
 3.141(452472285462...),z5 = 3.1415926(64850249...)

There are a few quadratic algorithms converging to the constant p. Quadratic convergence means that the number of correct digits is doubled at each iteration. The existence of such extraordinarily rapid algorithms using only basic operations (multiplication, addition, division) and square roots extraction is a remarkable fact.

#### 1.2.1  First algorithm

The first known quadratic algorithm was published in 1976 by Eugene Salamin in [6] and it was also discovered independently by Richard Brent the same year [5]. It uses the real AGM iteration.

Algorithm

 xk+1
 =
 12 ( xk+yk)
 yk+1
 =
 Ö xkyk
 ak+1
 =
 ak-2k+1( xk+12-yk+12)

with the initial values

 x0 = 1,y0 = 1Ö2 ,a0 = 12

then

 lim k® ¥ pk = lim k® ¥ 2xk2ak = p.

The first four iterations are

 p1
 =
 3.1(876...)
 p2
 =
 3.141(680...)
 p3
 =
 3.141592653(895...)
 p4
 =
 3.14159265358979323846(636...)

and computing p30 will provide more than a billion digits of p ! So the computation of only 30 square roots and less than a hundred multiplications is required to find a huge number of digits.

#### 1.2.2  Second algorithm

This second algorithm is due to Jonathan and Peter Borwein and was published in 1984 [1].

Algorithm

 xk+1
 =
1
2
æ
ç
ç
ç
è

Ö

xk

+ 1
 Ö xk
ö
÷
÷
÷
ø
 yk+1
 =
 Ö xk æç è yk+1yk+xk ö÷ ø
 ak+1
 =
 akyk+1 æç è xk+1+1yk+1+1 ö÷ ø

starting with

 x0 = Ö2,y0 = 0,a0 = 2+Ö2

then

 | ak-p| < 10-2k.

The rate of convergence is also quadratic and the first iterations are

 a1
 =
 3.14(260...)
 a2
 =
 3.1415926(609...)
 a3
 =
 3.141592653589793238(645...)

#### 1.2.3  Third algorithm

In [2], the following algorithm is given :

Algorithm

 yk+1
 =
 1- Ö 1-yk2

 1+ Ö 1-yk2
 ak+1
 =
 (1+yk+1)2ak-2k+1yk+1,

starting with

 y0 = 1/Ö2,a0 = 1/2,

then

 0 < ak- 1p < 16.2ke-p.2k.

The rate of convergence is again quadratic (the number of digits is doubled at each iteration). The first iterates are given by

 1/a1
 =
 2.9142135623730950488016887...
 1/a2
 =
 3.14(057...)
 1/a3
 =
 3.1415926(462...)
 1/a4
 =
 3.141592653589793238(279...)

and 1/a5 has 40 correct digits, 1/a6 has 83 correct digits, ...

### 1.3  Cubic

In [3], a cubic algorithm is described :

Algorithm

 xk+1
 =
 31+2(1-yk3)1/3
 yk+1
 =
 xk+1-12
 ak+1
 =
 xk+12ak-3k(xk+12-1),

with the initial values

 y0 = (Ö3-1)/2,a0 = 1/3,

then

 lim k® ¥ 1ak = p

and the rate of convergence is now cubic (the number of digits is tripled at each iteration). Again let's take a look to the first iterations

 1/a1
 =
 3.14159(058...)
 1/a2
 =
 3.141592653589793238462(359...)
and 1/a3 has 70 correct digits...

### 1.4  Quartic

Two iterations of the third quadratic algorithm gives a new iteration [2] :

Algorithm

 yk+1
 =
 1-( 1-yk4) 1/41+(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

 0 < ak- 1p < 16.4ke-2p.4k

and it has the impressive quartic rate of convergence (each successive iteration quadruples the number of correct digits). The first two iterates are

 1/a1
 =
 3.1415926(462...)
 1/a2
 =
 3.1415926535897932384626433832795028841971(146...)

This algorithm is among the most efficient and is often used with a quadratic one to compute billions of digits of the constant p (one is used for the computation and the other is used to check the result of the first one).

### 1.5  Quintic

In [4], this remarkable algorithm is given :

Algorithm

 yk+1
 =
 25yk(c+a/c+1)2
 a
 =
 5yk -1,b = (a-1)2+7,c = æç è 12 a(b+ Ö b2-4a3 ö÷ ø 1/5
 ak+1
 =
 yk2ak-5k æç è 12 (yk2-5) + Ö yk(yk2-2yk+5) ö÷ ø ,

starting with

 y0 = 5(Ö5-2),a0 = 1/2,

we have

 0 < ak- 1p < 16.5ke-p.5k

and the rate of convergence is quintic (the number of digits is multiplied by five at each iteration).

 1/a1
 =
 3.1415(369...)
 1/a2
 =
 3.141592653589793238462643383279(351...)

More complex algorithms with higher order are available like septic or nonic algorithms (in the nonic iteration the number of correct digits is multiplied by 9 at each step !).

### 1.6  Other algorithms

Starting with any approximation of p, the following algorithm was given by Beeler in 1972 and has a cubic convergence

 ak+1 = ak+sin( ak)

and for example using Archimedes' approximation

 a0 = 227

it gives

 a1
 =
 3.141592653(926...)
 a2
 =
 3.1415926535897932384626433832(858...)

and a3 has 88 correct digits,...

It's possible to find algorithms of this kind at any order and the following iterations have respectively quintic, septic, nonic, ... order of convergence

 ak+1
 =
 ak+sin( ak) + 16 sin3( ak)
 ak+1
 =
 ak+sin( ak) + 16 sin3( ak) + 340 sin5( ak)
 ak+1
 =
 ak+sin( ak) + 16 sin3( ak) + 340 sin5( ak) + 5112 sin7( ak)

the coefficients are those of the series expansion of the arcsine function. If one prefers to evaluate tangent functions we also have

 ak+1
 =
 ak-tan( ak)
 ak+1
 =
 ak-tan( ak) + 13 tan3( ak)
 ak+1
 =
 ak-tan( ak) + 13 tan3( ak) - 15 tan5( ak)
 ak+1
 =
 ak-tan( ak) + 13 tan3( ak) - 15 tan5( ak)+ 17 tan7( ak)

with respectively cubic, quintic, septic and nonic convergence and the coefficients are those of the expansion of the arctan function. For example, starting with the famous approximation a0 = 355/113, just one iteration of the nonic iterations gives 60 correct digits of p !. Similar formulae exist with cosines.

The major drawback of those algorithms is to evaluate sine or tangent function at a given accuracy which is not so easy.

## References

[1]
J.M. Borwein and P.B. Borwein, The Arithmetic-Geometric Mean and Fast Computation of Elementary Functions, SIAM review, (1984), vol. 26, p. 351-366

[2]
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)

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

[4]
J.M. Borwein, P.B. Borwein and D.H. Bailey, Ramanujan, Modular Equations, and Approximations to Pi or How to Compute One Billion Digits of Pi, The American Mathematical Monthly, (1989), vol. 96, p. 201-219

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

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

Back to The constant pi

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