**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-20^{th} 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 19^{th} 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 20^{th} CENTURY

__ 2.1. From ancient times till the 16^{th} century:__ The earliest values of p
including the "Biblical" value of 3, were almost certainly found by measurements. In the Egyptian

__ 2.2. From the 16^{th} till the 20^{th} century:__ The introduction of

__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 20^{th} CENTURY

In the beginning of the 20^{th} 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 18^{th} 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 a_{0}=1; b_{0}=1/Ö
2 and s_{0}=0.5 and calculate the iteration for k=1,2,...:

(8)

Then *p _{k} *converges to 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 a_{0} = 1/3 and s_{0}=(Ö
3-1)/2 and iterate:

(9)

then the series of *1/a _{k}* converges

1 term (k=1): r_{1} »
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 a_{0} = 6-4Ö
2 and y_{0} = Ö
2-1 and iterate [9]:

(10)

and *1/a _{k}* converges to p

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=n ^{2}* operations. With scientific terms: the bit complexity of multiplication is O(

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:=(x _{0}, x_{1}, x_{2}, ..., x_{n-1})* and

(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

(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=2 ^{m}*, then the discrete FFT can be evaluated in only

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 N^{th} 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.