PI, FOURIER TRANSFORM AND LUDOLPH VAN CEULEN

M.Vajta

Department of Mathematical Sciences

University of Twente

P.O.Box 217, 7500 AE Enschede

The Netherlands

e-mail: m.vajta@math.utwente.nl

 

 

ABSTRACT

The paper describes an interesting (and unexpected) application of the Fast Fourier transform in number theory. Calculating more and more decimals of p (first by hand and then from the mid-20th century, by digital computers) not only fascinated mathematicians from ancient times but kept them busy as well. They invented and applied hundreds of methods in the process but the known number of decimals remained only a couple of hundred as of the late 19th century. All that changed with the advent of the digital computers. And although digital computers made possible to calculate thousands of decimals, the underlying methods hardly changed and their convergence remained slow (linear). Until the 1970's. Then, in 1976, an innovative quadratic convergent formula (based on the method of algebraic-geometric mean) for the calculation of p was published independently by Brent [10] and Salamin [14]. After their breakthrough, the Borwein brothers soon developed cubically and quartically convergent algorithms [8,9]. In spite of the incredible fast convergence of these algorithms, it was the application of the Fast Fourier transform (for multiplication) which enhanced their efficiency and reduced computer time [2,12,15].

The author would like to dedicate this paper to the memory of Ludolph van Ceulen (1540-1610), who spent almost his whole life to calculate the first 35 decimals of p .

Keywords: approximation theory, Fast Fourier transform, elliptic function.

 

1. INTRODUCTION

On july 5, 2000 a very special ceremony took place in the St.Pieterskerk (St.Peter's Church) at Leiden, the Netherlands. A replica of the original tombstone of Ludolph van Ceulen was placed into the Church since the original disappeared [6, pp.51,13]. Ludolph van Ceulen (born 28 January 1540 in Hildesheim and died 31 December 1610 in Leiden) was a mathematician and fencing teacher. In 1600 Prince Maurits appointed him as one of the first hoogleraar wiskunde (professor of mathematics) at the University of Leiden. He dedicated almost all of his life to calculate more and more decimals of p . He published the 20 decimals in his book: Van de Circel in 1596. But he went on to calculate the first 35 decimals. As legend has it, the 35 decimals were engraved on his tombstone which later disappeared.

Now, 400 years later, the known number of digits of p exceeds 200 billion and increasing. It is not the number crunching, which is interesting, but the new methods developed in the last decades. Although the best algorithms seem to be simple in form but their computer implementation takes some time and ingenuity. The key element is an efficient multiplication method. And here we can find an unexpected application from signal processing: the Fast Fourier transform (FFT).

2. SHORT HISTORY TILL THE 20th CENTURY

2.1. From ancient times till the 16th century: The earliest values of p including the "Biblical" value of 3, were almost certainly found by measurements. In the Egyptian Rhind Papyrus (about 1650 BC) there is good evidence for p » 4(8/9)2 =3,16. The first theoretical calculation seems to have been carried out by Archimedes of Syracuse (287-212 BC). He used inscribed and circumscribed polygons. Applying a polygon with 96 side he obtained the approximation 310/71 < p < 31/7 (3,1408<p <3,1429). Archimedes did not have the advantage of an algebraic and trigonometric notation so he had to rely on pure geometrical means. Almost all of the methods from this period are based on the approximations of inscribed and circumscribed polygons. The highlight of this period was achieved by Ludolph van Ceulen (1540-1610) who determined the first 20 digits about 1596 based on a polygon with 60 x 229 (=32212254720) sides [7, pp.34]. Some years later he succeeded to establish the first 35 digits [6, pp.102,13]. The Germans were so much impressed by van Ceulen's achievement that they began to call p as the Ludolph's number.

2.2. From the 16th till the 20th century: The introduction of calculus in the 16th century made a large number of new algorithms in the form of infinite series, products or continued fractions possible. François Viète, a lawyer and amateur mathematician, was the first who expressed p by an infinite product in 1593. John Wallis also developed an infinite product but in much simpler form. Newton used his own formulae to calculate the first 15 decimals of p . When Gregory discovered the arctangent series in 1671 it led to a number of new algorithms [6:pp.92]. In 1706 Machin used Gregory's series and his formula (named after him) to calculate 100 decimals of p [2]. Euler invented dozens of new formulae for p [6:pp.112]. Without going into details we simply present some of the most known formulae from this period:

F.Viète (ca 1593): (1)

J.Wallis (ca 1655): (2)

I.Newton (ca 1666): (3)

Gregory (1671): (4)

Machine (1706): (5)

L.Euler (ca 1736): (6)

As for the nature of p , Lambert proved it to be irrational in 1761. It took more then 100 years till Lindemann finally proved in 1882 that p is in fact transcendental putting the problem of squaring the circle to rest.

 

 

3. THE 20th CENTURY

In the beginning of the 20th century, there appeared some new and very interesting theoretical result by Ramanujan (based on elliptic and modular functions) [9]. Besides the theoretical development it was the introduction of digital computers which made possible to calculate an unbelievable large number of decimals of p . In the first 4000 years mathematicians could only calculate (by hand) the first 707 digits. As a test of the first digital computer (ENIAC) John von Neumann proposed to calculate the value of p . The first calculation by ENIAC in 1947 already provided 2000 digits. The team of scientists needed 70 hours computer time including programming [6, pp.277].

With the computer a new era began. As the speed and memory of computers grew, so increased the known number of digits. In 1961 D.Shanks and J.W.Wrench,Jr. break the 100000 decimals barrier on an IBM 7090 computer in less then 9 hours. But the formula used was still a variant of the arctan formula due to Störmer [6, pp.576]:

(7)

Interestingly, almost all calculations from the beginning of the 18th century until the early 1970's have relied on one or an other form of the Machin's (arctan) formula. Since the 60' the increase in number of decimals has been quite remarkable. Although Shanks and Wrench predicted in 1961 that to calculate the first 1 million digit should last for month (in computer time), nowadays it takes only minutes to calculate a couple of millions decimals on a fast Pentium PC. In 1996 the Chudnovsky brothers crossed the 1 billion limit. How could they do that? With faster and faster new algorithms and with an efficient multiplication method by the Fast Fourier transform.

 

4.1. The Brent-Salamin algorithm:

One of the most crucial drawbacks of the known algorithms is their slow (linear) convergence. It was therefore a milestone in the history of p research when Eugene Salamin and Richard Brent developed (independently from each other) a dramatically new algorithm in the 70's [10,14]. Their algorithm is based on the arithmetic-geometric mean (AGM) known already to Gauss [8]. The arithmetic-geometric mean was the basis of Gauss' method for the calculation of elliptic integrals. With the help of the elliptic integral relation of Legendre, p can be expressed in terms of the arithmetic-geometric mean and the resulting algorithm is quadratically fast [12]. The Brent-Salamin algorithm is the following: set a0=1; b0=1/Ö 2 and s0=0.5 and calculate the iteration for k=1,2,...:

(8)

Then pk converges to p quadratically. The first 5 terms of the Brent-Salamin algorithm give 1,3,9,20 and 42 decimals of p :

1 term (k=1): p » 3,1 ...

2 terms (k=2): p » 3,141 ...

3 terms (k=3): p » 3,141592653 ...

4 terms (k=4): p » 3,1415926535 8979323846 ...

5 terms (k=5): p » 3,1415926535 8979323846 2643383279 5028841971 69 ...

Quadratic converges means that with each new term the number of digits doubles! For a while it seemed there was no way to develop faster algorithms. But not for long. The results of Brent and Salamin gave a new impetus to the p research. Built on the same body of mathematics Jonathan and Peter Borwein introduced an even faster algorithm in 1985.

 

4.2. Borweins' cubically convergent algorithm:

Jonathan and Peter Borwein developed an even faster algorithms to approximate p in the 1980's. They also applied the arithmetic-geometric mean (AGM) and transformation theory of elliptic integral and modular equations [8,9]. Their cubically convergent algorithm is as follows: set a0 = 1/3 and s0=(Ö 3-1)/2 and iterate:

(9)

then the series of 1/ak converges cubically to p . Thus provides the first three terms 5,21 and 57 digits:

1 term (k=1): r1 » 3.14159 ...

2 terms (k=2): p » 3.1415926535 8979323846 2 ...

3 terms (k=3): p » 3.1415926535 8979323846 2643383279 5028841971

6939937510 5820974 ...

 

4.3. Borweins' quarticall convergent algorithm:

It was a great achievement of the Borwein's to push further and establish a quartically convergent algorithm. Set a0 = 6-4Ö 2 and y0 = Ö 2-1 and iterate [9]:

(10)

and 1/ak converges to p quartically. The first terms provides:

1 term (k=1): p » 3.1415926 ...

2 terms (k=2): p » 3.1415926535 8979323846 2643383279 5028841971

2 terms (k=2): p » 3.1415926535 8979323846 2643383279 5028841971

6939937510 5820974944 5923078164 0628620899

8628034825 3421170679 8214808651 3282306647

0938446095 ...

The first 2 terms already gives 40 decimals, the first three terms gives 170 decimals! Only if Ludolph van Ceulen had known this method! In fact, this series converge so fast that taking only the first 15 terms provides 2 billions digits of p ! Bailey applied this algorithm in his record-breaking calculation in 1987 [2].

It has been a long way to come from Machin's formulae to the quartically convergent Borwein-Borwein algorithm. But to increase the efficiency of the calculations one need better (i.e. faster) numerical methods as well. It seems simple to program algorithm (9) but how can we multiply two large numbers efficiently?

 

5. THE KEY INGREDIENT: FAST FOURIER TRANSFORM

Algorithm (8), (9) or (10) does not seem difficult to be implemented on a digital computer. We must note, however, that all operations must be correct up to the required number of digits plus m (30 as the guard digits for millions of decimal digit calculations) [12]. So we can see, that an efficient multiplication method is a key element in all algorithms. We learned how to multiply two numbers already in the elementary school and we know that to multiply two n-digit numbers we need nxn=n2 operations. With scientific terms: the bit complexity of multiplication is O(n2). It seems all pretty simple but when we have to multiply two numbers with millions of digits, then it is another story. So we may ask: isn't there any better way to multiply large (i.e. long) numbers? The unexpected answer came in 1971. Then, Schönhage and Strassen showed that it is possible to multiply two n-digit integers with bit complexity O(nlogn/loglogn) [9,15].

Their method based on the application of the Fourier transform. Fourier transform? Fourier transform has long been known to mathematicians and engineers alike. Its application ranges from heat equations to sound engineering to speech recognition. But in number theory? How could that be? Here is a short answer. Suppose we have to multiply two n-digit numbers x:=(x0, x1, x2, ..., xn-1) and y:=(y0, y1, y2, ..., yn-1). Then a key observation may be made: the product sequence z:=(z0, z1, z2, ..., z2n-1) of x and y is precisely the discrete convolution C(x,y):

(11)

where the subscript k-j is to be interpreted as k-j+N if k-j is negative. Now we apply the well known discrete Fourier Transform [11]. First, extend x and y to length N=2n by appending zeros at the end of each. Let F(x) denote the discrete Fourier transform of the sequence x, and let F-1(x) denote the inverse discrete Fourier transform of x [2,9,11,12]:

(12)

Then the convolution theorem states, that the Fourier transform of a convolution product is the ordinary product of the Fourier transforms:

(13)

Thus the enirely multiplication pyramid z can be obtained by performing two forward discrete Fourier transforms, one vector complex multiplication and one inverse transform, each of length N=2n.

One must realize that it is the discrete Fast Fourier transform (FFT) which makes this scheme work. In particular, if N=2m, then the discrete FFT can be evaluated in only 5m2m arithmetic operations.

There are of course several "tricks" in implementing FFT based multiplication. One usual trick is to utilize the fact that the input data vectors x and y and the result vector z are purely real.

One other variation relies on the fact that the FFT can be applied in any number field in which there exists a primitive Nth root of unity. This requirements holds for the field of integers modulo p, where p is a prime of the form p=kN+1 [9]. The advantage to use a prime modulus field (instead of the field of complex numbers) is that there are no round-off errors (since all computations are exact). Some further details concerning the implementations can be found in [2,9,12].

Calculation of Reciprocals and Square Roots Do we need something else besides a very fast multiplication method? Not really. For we can determine the inverse of a number 1/y or its square root Ö y or 1/Ö y by the Newton's quadratically convergent method:

1/y: →

√y: → (14)

1/√y: →

By applying the FFT method to multiply large numbers, it is possible to accelerate the computations dramatically. In fact, in almost all record-breaking high-performance multi-precision computer programs recently, some variation of the FFT method has been applied [2,9,12]. The current record of number of digits is more than 200 billion. So we can conclude, that the successful combination of number theory and computer science (via the FFT) made these records possible.

 

 

CONCLUSIONS

Some methods of approximating the value of p are presented. With new methods developed in the last decades it became possible to calculate billions of decimals of p . But to program the fastest algorithms of the Borweins' one must have efficient multiplication methods as well. It is the discrete Fast Fourier transform, which made fast multiplication of very long numbers possible. Almost all the current records apply one or another version of FFT multiplication. In recent years, the computation of the expansion of p has assumed the role as a standard test of computer integrity. Finally, it is interesting to note that the p research enjoys considerable wide attention: popular books appeared on the subject [5,7] and even a new cologne named Pi is marketed by Givenchy.

 

REFERENCES

[1] ABRAMOWITZ,M. - I.A.STEGUN (Eds.): Handbook of Mathematical Functions,

Dover Publications, Inc., New York, 1972.

[2] BAILEY,D.H: The Computation of p to 29.360.000 Decimal Digits Using Borweins'

Quartically Convergent Algorithm,

Mathematics of Computation, 50, No.181, (1988), pp.283-296.

[3] BAILEY,D.H., J.BORWEIN, P.BORWEIN and S.PLOUFFE: The Quest for PI,

Mathematical Intelligencer, 19, No.1, (1997), pp.50-57.

[4] BAILEY,D.H., P.BORWEIN and S.PLOUFFE: On the Rapid Computation of Various

Polylogarithmic Constants,

in PI: A Source Book, Springer-Verlag, New York, 1997. pp.663-676.

[5] BECKMAN,P: A history of p (pi),

Golem Press (St.Martin's Press), New York, 1971.

[6] BERGGREN,L., J.BORWEIN and P.BORWEIN: PI: A Source Book,

Springer-Verlag, New York, 1997.

[7] BLATNER,D: The Joy of p ,

Penguin Books, London, 1997.

[8] BORWEIN,J.M. and P.B.BORWEIN: Pi and the AGM - A Study in Analytic Number Theory

and Computational Complexity,

John Wiley, New York, 1987.

[9] BORWEIN,J.M., 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, 96, No.3, (1989), pp.201-219.

[10] BRENT,R.P: Fast Multiple-Precision Evaluation of Elementary Functions,

Journal of the Association of Computing Machinery, 23, No.2, (1976), pp.242-251.

[11] KAMEN,E.W. and B.S.HECK: Fundamentals of Signals and Systems using Matlab,

Prentice-Hall, Inc., Upper Saddle River, New Jersey, 1997.

[12] KANADA,Y.: Vectorization of Multiple-Precision Arithmetic Program and 201.326.000 Decimal

Digits of p Calculation,

Scientific American, 1988, also in PI: A Source Book, pp.576.

[13] "PI in de Pieterskerk" (Pi in the ST.Peter's Church),

Nieuw Archief voor Wiskunde (NAW), 5e serie, No.2, (juni 2000), pp.114-115.

[14] SALAMIN,E: Computation of p Using Arithmetic-Geometric Mean,

Mathematics of Computation, 30, No.135, (1976), pp.565-570.

[15] SCHÖNHAGE,A. and V.STRASSEN: Schnelle Multiplikation grosser Zahlen,

Computing (Arch. Elektr. Rechnen), 7, (1971), pp.281-292.