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
=   
 


1
2
 
  
 
 


1
2
+ 1
2
  
 


1
2
 
 
  
 
 


1
2
+ 1
2
  
 


1
2
+ 1
2


 

[1/2]
 
 
 
...

an algorithm may be written as

xk+1
=
xkyk
yk+1
=
  
 


2yk
yk+1
 

and

x0
=
1
y0
=
2

then


lim
k  
xk = p
2
.

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
=
2xkyk
xk+yk
yk+1
=


 

xk+1yk
 

starting with

x0 = 23,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  


1
3
xk+ 2
3
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...)

1.2  Quadratic

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
=
1
2
( xk+yk)
yk+1
=


 

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

with the initial values

x0 = 1,y0 = 1
2
,a0 = 1
2

then


lim
k  
pk =
lim
k  
2xk2
ak
= 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+1
yk+xk


ak+1
=
akyk+1

xk+1+1
yk+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- 1
p
< 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
=
3
1+2(1-yk3)1/3
yk+1
=
xk+1-1
2
ak+1
=
xk+12ak-3k(xk+12-1),

with the initial values

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

then


lim
k  
1
ak
= 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/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-42,

then

0 < ak- 1
p
< 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
=
25
yk(c+a/c+1)2
a
=
5
yk
-1,b = (a-1)2+7,c =

1
2
a(b+

 

b2-4a3
 


1/5

 
ak+1
=
yk2ak-5k

1
2
(yk2-5) +

 

yk(yk2-2yk+5)
 


,

starting with

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

we have

0 < ak- 1
p
< 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 = 22
7

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) + 1
6
sin3( ak)
ak+1
=
ak+sin( ak) + 1
6
sin3( ak) + 3
40
sin5( ak)
ak+1
=
ak+sin( ak) + 1
6
sin3( ak) + 3
40
sin5( ak) + 5
112
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) + 1
3
tan3( ak)
ak+1
=
ak-tan( ak) + 1
3
tan3( ak) - 1
5
tan5( ak)
ak+1
=
ak-tan( ak) + 1
3
tan3( ak) - 1
5
tan5( ak)+ 1
7
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.