Ramanujan, Modular Equations, and Approximations to Pi or How to compute One Billion Digits of Pi

The year 1987 was the centenary of Ramanujan’s birth. He died in 1920 Had he not died so young, his presence in modern mathematics might be more immediately felt. Had he lived to have access to powerful algebraic manipulation software. such as MACSYMA, who knows how much more spectacular his already astonishing career might have been.


1.
Introduction. There is a close and beautiful connection between the transformation theory for elliptic integrals and the very rapid approximation of pi. This connection was first made explicit by Ramanujan in his 1914 paper "Modular Equations and Approximations to 7r" [26]. We might emphasize that Algorithms 1 and 2 are not to be found in Ramanujan's work, indeed no recursive approximation of g" is considered, but as we shall see they are intimately related to his analysis. and an converges to 1/ST quartically (that is, with order four).
Algorithm 2. Let so = 5(V5 -2) and a0 -1/2. Then 0 < --< 16 5 ne-5n., 7T and an converges to 1/ST quintically (that is, with order five). Each additional term in Sum 1 adds roughly eight digits, each additional iteration of Algorithm 1 quadruples the number of correct digits, while each additional iteration of Algorithm 2 quintuples the number of correct digits. Thus a mere thirteen iterations of Algorithm 2 provide in excess of one billion decimal digits of pi. In general, for us, pth-order convergence of a sequence {tan) to a means that an tends to a and that laxn+1-al < Claxn -aP for some constant C > 0. Algorithm 1 is arguably the most efficient algorithm currently known for the extended precision calculation of pi. While the rates of convergence are impressive, it is the subtle and thoroughly nontransparent nature of these results and the beauty of the underlying mathematics that intrigue us most.
Watson [37], commenting on certain formulae of Ramanuj an, talks of a thrill which is indistinguishable from the thrill which I feel when I enter the Sagrestia Nuovo of the Capella Medici and see before me the austere beauty of the four statues representing "Day," "Night," "Evening," and "Dawn" which Michelangelo has set over the tomb of Giuliano de'Medici and Lorenzo de 'Medici.
Sum 1 is directly due to Ramanujan and appears in [26]. It rests on a modular identity of order 58 and, like much of Ramanujan's work, appears without proof and with only scanty motivation. The first complete derivation we know of appears in [11]. Algorithms 1 and 2 are based on modular identities of orders 4 and 5, respectively. The underlying quintic modular identity in Algorithm 2 (the relation for sn) is also due to Ramanujan, though the first proof is due to Berndt and will appear in [7]. One intention in writing this article is to explain the genesis of Sum 1 and of Algorithms 1 and 2. It is not possible to give a short self-contained account without assuming an unusual degree of familiarity with modular function theory. Also, parts of the derivation involve considerable algebraic calculation and may most easily be done with the aid of a symbol manipulation package (MACSYMA, MAPLE, REDUCE, etc.). We hope however to give a taste of methods involved. The full details are available in [11].
A second intention is very briefly to describe the role of these and related approximations in the recent extended precision calculations of pi. In part this entails a short discussion of the complexity and implementation of such calculations. This centers on a discussion of multiplication by fast Fourier transform methods. Of considerable related interest is the fact that these algorithms for ST are provably close to the theoretical optimum.
2. The State of Our Current Ignorance. Pi is almost certainly the most natural of the transcendental numbers, arising as the circumference of a circle of unit diameter. Thus, it is not surprising that its properties have been studied for some twenty-five hundred years. What is surprising is how little we actually know.
We know that ST is irrational, and have known this since Lambert's proof of 1771 (see [5]). We have known that g is transcendental since Lindemann's proof of 1882 [23]. We also know that 7T  We know that e ' is transcendental. This follows by noting that es = (1 1)-i and applying the Gelfond-Schneider theorem [4]. We know that ST + log 2 + V2log 3 is transcendental. This result is a consequence of the work that won Baker a Fields Medal in 1970. And we know a few more than the first two hundred million digits of the decimal expansion for 7T (Kanada, see Section 3). The state of our ignorance is more profound. We do not know whether such basic constants as ST + e, 7r/e, or log ST are irrational, let alone transcendental. The best we can say about these three particular constants is that they cannot satisfy any polynomial of degree eight or less with integer coefficients of average size less than 109 [3]. This is a consequence of some recent computations employing the [March Ferguson-Forcade algorithm [17]. We don't know anything of consequence about the single continued fraction of pi, except (numerically) the first 17 million terms, which Gosper computed in 1985 using Sum 1. Likewise, apart from listing the first many millions of digits of Tr, we know virtually nothing about the decimal expansion of 7T. It is possible, albeit not a good bet, that all but finitely many of the decimal digits of pi are in fact O's and l's. Carl Sagan's recent novel Contact rests on a similar possibility. Questions concerning the normality of or the distribution of digits of particular transcendentals such as X appear completely beyond the scope of current mathematical techniques. The evidence from analysis of the first thirty million digits is that they are very uniformly distributed [2]. The next one hundred and seventy million digits apparently contain no surprises.
In part we perhaps settle for computing digits of 7r because there is little else we can currently do. We would be amiss, however, if we did not emphasize that the extended precision calculation of pi has substantial application as a test of the "global integrity" of a supercomputer. The extended precision calculations described in Section 3 uncovered hardware errors which had to be corrected before those calculations could be successfully run. Such calculations, implemented as in Section 4, are apparently now used routinely to check supercomputers before they leave the factory. A large-scale calculation of pi is entirely unforgiving; it soaks into all parts of the machine and a single bit awry leaves detectable consequences.

Matters Computational
I am ashamed to tell you to how many figures I carried these calculations, having no other business at the time.
Isaac Newton Newton's embarrassment at having computed 15 digits, which he did using the arcsinlike formula -x2 dx, is indicative both of the spirit in which people calculate digits and the fact that a surprising number of people have succumbed to the temptation [5]. The history of efforts to determine an accurate value for the constant we now know as X is almost as long as the history of civilization itself. By 2000 B.C. both the Babylonians and the Egyptians knew ST to nearly two decimal places. The Babylonians used, among others, the value 3 1/8 and the Egyptians used 3 13/81. Not all ancient societies were as accurate, however-nearly 1500 years later the Hebrews were perhaps still content to use the value 3, as the following quote suggests.
Also, he made a molten sea of ten cubits from brim to brim, round in compass, and five cubits the height thereof; and a line of thirty cubits did compass it round about.
Old While hardly better than estimates one could get with a ruler, this is the first method that can be used to generate an arbitrary number of digits, and to a nonnumerical mathematician perhaps the problem ends here. Variations on this theme provided the basis for virtually all calculations of ST  This particular formula is due to Machin and was employed by him to compute 100 digits of g in 1706. Variations on this second theme are the basis of all the calculations done until the 1970's including William Shanks' monumental handcalculation of 527 digits. In the introduction to his book [32], which presents this calculation, Shanks writes: Towards the close of the year 1850 the Author first formed the design of rectifying the circle to upwards of 300 places of decimals. He was fully aware at that time, that the accomplishment of his purpose would add little or nothing to his fame as a Mathematician though it might as a Computer; nor would it be productive of anything in the shape of pecuniary recompense.
Shanks actually attempted to hand-calculate 707 digits but a mistake crept in at the 527th digit. This went unnoticed until 1945, when D. Ferguson, in one of the last "nondigital" calculations, computed 530 digits. Even with machine calculations mistakes occur, so most record-setting calculations are done twice-by sufficiently different methods.
The advent of computers has greatly increased the scope and decreased the toil of such calculations. Metropolis, Reitwieser, and von Neumann computed and analyzed 2037 digits using Machin's formula on ENIAC in 1949. In 1961, Dan Shanks and Wrench calculated 100,000 digits on an IBM 7090 [31]. By 1973, still using Machin-like arctan expansions, the million digit mark was passed by Guillard and Bouyer on a CDC 7600.
[March iii) The third technique, based on the transformation theory of elliptic integrals, provides the algorithms for the most recent set of computations. The most recent records are due separately to Gosper, Bailey, and Kanada. Gosper in 1985 calculated over 17 million digits (in fact over 17 million terms of the continued fraction) using a carefully orchestrated evaluation of Sum 1.
Bailey in January 1986 computed over 29 million digits using Algorithm 1 on a Cray 2 [2]. Kanada, using a related quadratic algorithm (due in basis to Gauss and made explicit by Brent [12] and Salamin [27]) and using Algorithm 1 for a check, verified 33,554,000 digits. This employed a HITACHI S-810/20, took roughly eight hours, and was completed in September of 1986. In January 1987 Kanada extended his computation to 227 decimal places of r and the hundred million digit mark had been passed. The calculation took roughly a day and a half on a NEC SX2 machine. Kanada's most recent feat (Jan. 1988) was to compute 201,326,000 digits, which required only six hours on a new Hitachi S-820 supercomputer. Within the next few years many hundreds of millions of digits will no doubt have been similarly computed. Further discussion of the history of the computation of pi may be found in [5] and [9].
4. Complexity Concerns. One of the interesting morals from theoretical computer science is that many familiar algorithms are far from optimal. In order to be more precise we introduce the notion of bit complexity. Bit complexity counts the number of single operations required to complete an algorithm. The single-digit operations we count are +, -, X. (We could, if we wished, introduce storage and logical comparison into the count. This, however, doesn't affect the order of growth of the algorithms in which we are interested.) This is a good measure of time on a serial machine. Thus, addition of two n-digit integers by the usual method has bit complexity 0(n), and straightforward uniqueness considerations show this to be asymptotically best possible.
Multiplication is a different story. Usual multiplication of two n-digit integers has bit complexity 0(n2) and no better. However, it is possible to multiply two n-digit integers with complexity 0(n (log n)(log log n)). This result is due to Schonhage and Strassen and dates from 1971 [29]. It provides the best bound known for multiplication. No multiplication can have speed better than 0(n). Unhappily, more exact results aren't available.
The original observation that a faster than 0(n2) multiplication is possible was due to Karatsuba in 1962. Observe that (a + blOn)(c + dlon) = ac + [(ab)(cd) -ac -bdIlOn + bd102n, and thus multiplication of two 2n-digit integers can be reduced to three multiplications of n-digit integers and a few extra additions. (Of course multiplication by 10" is just a shift of the decimal point.) If one now proceeds recursively one produces a multiplication with bit complexity 0(n log2 3) Note that 1og23 = 1.58... < 2.
We denote by M(n) the bit complexity of multiplying two n-digit integers together by any method that is at least as fast as usual multiplication.
The trick to implementing high precision arithmetic is to get the multiplication right. Division and root extraction piggyback off multiplication using Newton's method. One may use the iteration Xk+1 = 2Xk ky to compute l/y and the iteration Xk?1 =(Xk + ) 2 Xk to compute Vyi. One may also compute 1/ jy from xk(3 -yX2) Xk+1' 2 and so avoid divisions in the computation of Jy. Not only do these iterations converge quadratically but, because Newton's method is self-correcting (a slight perturbation in Xk does not change the limit), it is possible at the kth stage to work only to precision 2k. If division and root extraction are so implemented, they both have bit complexity O(M(n)), in the sense that n-digit input produces n-digit accuracy in a time bounded by a constant times the speed of multiplication. This extends in the obvious way to the solution of any algebraic equation, with the startling conclusion that every algebraic number can be computed (to n-digit accuracy) with bit complexity O(M(n)). Writing down n-digits of F/ or 3ff is (up to a constant) no more complicated than multiplication.
The Schonhage-Strassen multiplication is hard to implement. However, a multiplication with complexity O((log n)2+"n) based on an ordinary complex (floating point) fast Fourier transform is reasonably straightforward. This is Kanada's approach, and the recent records all rely critically on some variations of this technique.
To see how the fast Fourier transform may be used to accelerate multiplication, let x = (x0, xl, x2,..., xn-1) and y = (yo, Yi, Y2 .' . y,n -) be the representations of two high-precision numbers in some radix b. The radix b is usually selected to be some power of 2 or 10 whose square is less than the largest integer exactly representable as an ordinary floating-point number on the computer being used. Then, except for releasing each "carry," the product z := (z0, Zl, Z2,**, Z2n-1) of x and y may be written as zo = xoyo z1 = xoYl + xlyo Z2n-3 = Xn-lYn-2 + Xn-2Yn-I where the subscript kj is to be interpreted as kj + N if kj is negative. Now a well-known result of Fourier analysis may be applied. Let F(x) denote the discrete Fourier transform of the sequence x, and let F-(x) denote the inverse discrete Fourier transform of x: Then the "convolution theorem," whose proof is a straightforward exercise, states that or, expressed another way, Thus the entire 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. Once the real parts of the resulting complex numbers have been rounded to the nearest integer, the final multiprecision product may be obtained by merely releasing the carries modulo b. This may be done by starting at the end of the z vector and working backward, as in elementary school arithmetic, or by applying other schemes suitable for vector processing on more sophisticated computers. A straightforward implementation of the above procedure would not result in any computational savings-in fact, it would be several times more costly than the usual "schoolboy" scheme. The reason this scheme is used is that the discrete Fourier transform may be computed much more rapidly using some variation of the well-known "fast Fourier transform" (FFT) algorithm [13]. In particular, if N = 2', then the discrete Fourier transform may be evaluated in only 5m2m arithmetic operations using an FFT. Direct application of the definition of the discrete Fourier transform would require 22m+3 floating-point arithmetic operations, even if it is assumed that all powers of e-2,i/N have been precalculated. This is the basic scheme for high-speed multiprecision multiplication. Many details of efficient implementations have been omitted. For example, it is possible to take advantage of the fact that the input sequences x and y and the output sequence z are all purely real numbers, and thereby sharply reduce the operation count. Also, it is possible to dispense with complex numbers altogether in favor of performing computations in fields of integers modulo large prime numbers. Interested readers are referred to [2], [8], [13], and [22].
When the costs of all the constituent operations, using the best known techniques, are totalled both Algorithms 1 and 2 compute n digits of IT with bit complexity O( M(n )log n), and use O(log n) full precision operations.
The bit complexity for Sum 1, or for v using any of the arctan expansions, is between 0((log n)2M(n)) and 0(nM(n)) depending on implementation. In each case, one is required to sum 0(n) terms of the appropriate series. Done naively, one obtains the latter bound. If the calculation is carefully orchestrated so that the terms are grouped to grow evenly in size (as rational numbers) then one can achieve the former bound, but with no corresponding reduction in the number of operations.
The Archimedean iteration of section 2 converges like 1/4" so in excess of n iterations are needed for n-digit accuracy, and the bit complexity is 0(nM(n)).
Almost any familiar transcendental number such as e, y, '(3), or Catalan's constant (presuming the last three to be nonalgebraic) can be computed with bit complexity O((log n)M(n)) or O((log n)2M(n)). None of these numbers is known to be computable essentially any faster than this. In light of the previous observation that algebraic numbers are all computable with bit complexity O(M(n)), a proof that X cannot be computed with this speed would imply the transcendence of ST. It would, in fact, imply more, as there are transcendental numbers which have complexity O(M(n)). An example is 0.10100100001....
It is also reasonable to speculate that computing the nth digit of ST is not very much easier than computing all the first n digits. We think it very probable that computing the nth digit of ST cannot be 0(n).

The Miracle of Theta Functions
When I was a student, abelian functions were, as an effect of the Jacobian tradition, considered the uncontested summit of mathematics, and each of us was ambitious to make progress in this field. And now? The younger generation hardly knows abelian functions.
Felix Klein [21] Felix Klein's lament from a hundred years ago has an uncomfortable timelessness to it. Sadly, it is now possible never to see what Bochner referred to as " the miracle of the theta functions" in an entire university mathematics program. A small piece of this miracle is required here [6], [11], [28]. First some standard notations. The The first remarkable identity is Legendre's relation namely [March (for 0 < k < 1), which is pivotal in relating these quantities to pi. We also need to define two Jacobian theta functions 00 (2 (q) = E q(n+1/2)2 The definitions we have given above are not complete. We will be more precise in our discussion of X. One might bear in mind that much of the theory for X holds in considerably greater generality.
The fundamental region F we associate with X is the set of complex numbers What makes X a X-modular function is the fact that X is meromorphic in {im t > 0) and that X(w(t)) X= (t) for all w in the X-group, plus the fact that X tends to a definite limit (possibly infinite) as t tends to a vertex of the fundamental region (one of the three points (0, -1), (0,0 ), (i, oo)). Here we only allow convergence from within the fundamental region. Now some of the miracle of modular functions can be described. Largely because every point in the upper half plane is the image of a point in F under an element of the X-group, one can deduce that any X-modular function that is bounded on F is constant. Slightly further into the theory, but relying on the above, is the result that any two modular functions are algebraically related, and resting on this, but further again into the field, is the following remarkable result. Recall that q is given by (5.9). THEOREM 1. Let z be a primitive pth root of unity for p an odd prime. Consider the pth order modular equation for X as defined by Then the function Wp is a polynomial in x and X (independent of q), which has integer coefficients and is of degree p + 1 in both x and X.
The modular equation for X usually has a simpler form in the associated variables u = x1/8 and v -= X/8. In this form the 5th-order modular equation is given by and turns out to be a rational function of k(qP) and k(q).
The function X(t) is 1-1 on F and has a well-defined inverse, XA1, with branch points only at 0, 1 and x. This can be used to provide a one line proof of the " big" Picard theorem that a nonconstant entire function misses at most one value (as does exp). Indeed, suppose g is an entire function and that it is never zero or one; then exp(X-1(g(z))) is a bounded entire function and is hence constant. Littlewood suggested that, at the right point in history, the above would have been a strong candidate for a 'one line doctoral thesis'.
6. Ramanujan's Solvable Modular Equations. Hardy [19] commenting on Ramanujan's work on elliptic and modular functions says It is here that both the profundity and limitations of Ramanujan 's knowledge stand out most sharply.
We present only one of Ramanujan's modular equations.

)3 (q)
where for i = 1 and 2 ri = x(y y2 4x3) with 583(q 5) _1( ) X := and y := (x -1)2 + 7. This is a slightly rewritten form of entry 12(iii) of Chapter 19 of Ramanujan's Second Notebook (see [7], where Berndt's proofs may be studied). One can think of Ramanujan's quintic modular equation as an equation in the multiplier Mp of (5.13). The initial surprise is that it is solvable. The quintic modular relation for A, W5, and the related equation for A18, 25, are both nonsolvable. The Galois group of the sixth-degree equation Q (see (5.12)) over Q(v) is As and is nonsolvable. Indeed both Hermite and Kronecker showed, in the middle of the last century, that the solution of a general quintic may be effected in terms of the solution of the 5th-order modular equation (5.12) and the roots may thus be given in terms of the theta functions.
In fact, in general, the Galois group for Wp of (5.11) has order p( p + 1)( p -1) and is never solvable for p > 5. The group is quite easy to compute, it is generated by two permutations. If q := eizt then -T + 2 and T (2+ 1) are both elements of the A-group and induce permutations on the Ai of Theorem 1. For any fixed p, one can use the q-expansion of (5.10) to compute the effect of these transformations on the Ai, and can thus easily write down the Galois group. While Wp is not solvable over Q(A), it is solvable over Q(A, A0). Note that A0 is a root of Wp. It is of degree p + 1 because Wp is irreducible. Thus the Galois group for Wp over Q(A, A0) has order p(p -1). For p = 5, 7, and 11 this gives groups of order 20, 42, and 110, respectively, which are obviously solvable and, in fact, for general primes, the construction always produces a solvable group.
From (5.8) and (5.10) one sees that Ramanujan's modular equation can be rewritten to give A5 solvable in terms of A0 and A. Thus, we can hope to find an explicit solvable relation for Ap in terms of A and A0. For p = 3, Wp is of degree 4 and is, of course, solvable. For p = 7, Ramanujan again helps us out, by providing a solvable seventh-order modular identity for the closely related eta function defined by 00 (q) := qq2 11 (1 - The first interesting prime for which an explicit solvable form is not known is the "endecadic" (p = 11) case. We consider only prime values because for nonprime values the modular equation factors.
[March This leads to the interesting problem of mechanically constructing these equations. In principle, and to some extent in practice, this is a purely computational problem. Modular equations can be computed fairly easily from (5.11) and even more easily in the associated variables u and v. Because one knows a priori bounds on the size of the (integer) coefficients of the equations one can perform these calculations exactly. The coefficients of the equation, in the variables u and v, grow at most like 2n. (See [11].) Computing the solvable forms and the associated computational problems are a little more intricate-though still in principle entirely mechanical. A word of caution however: in the variables u and v the endecadic modular equation has largest coefficient 165, a three digit integer. The endecadic modular equation for the intimately related function J (Klein's absolute invariant) has coefficients as large as 27090964785531389931563200281035226311929052227303 x 29231952011253. It is, therefore, one thing to solve these equations, it is entirely another matter to present them with the economy of Ramanujan.
The paucity of Ramanujan's background in complex analysis and group theory leaves open to speculation Ramanujan's methods. The proofs given by Berndt are difficult. In the seventh-order case, Berndt was aided by MACSYMA-a sophisticated algebraic manipulation package. Berndt comments after giving the proof of various seventh-order modular identities: Of course, the proof that we have given is quite unsatisfactory because it is a verification that could not have been achieved without knowledge of the result. Ramanujan obviously possessed a more natural, transparent, and ingenious proof. (va(r)) K] (7.2) We have suppressed, and will continue to suppress, the k variable. With (5.6) and (5.7) at hand we can write a q-expansion for a, namely, and we can see that as r tends to infinity q = e-'IW tends to zero and a(r) tends to 1/7. In fact a(r) 1 ( 1) -e (7.4) The key now is iteratively to calculate a. This is the content of the next theorem. where represents the full derivative of Mp with respect to ko. In particular, a is algebraic for rational arguments.
We know that K(k1) is related via MP to K(k) and we know that E(k) is related via differentiation to K. (See (5.7) and (5.13).) Note that q -* qP corresponds to r -p2r. Thus from (7.2) some relation like that of the above theorem must exist. The actual derivation requires some careful algebraic manipulation.  Here (c)" is the rising factorial: (C) := c(c + 1)(c + 2) ... (c + n-1).
Some of the ingredients for the proof of Theorem 5, which are detailed in [11], are the following. Our first step is to write (7.2) as a sum after replacing the E by K and dK/dk using (5.7). One then uses an identity of Clausen's which allows one to write the square of a hypergeometric function 2F1 in terms of a generalized hypergeometric 3F2, namely, for all k one has   [26]; is presented by Watson in a long series of papers commencing with [36]; and some more recent derivations are given in [11] and [30]. An inspection of Theorem 5 shows that all the constants in Series 1 are determined from g58. Knowing a is equivalent to determining that the number 1103 is correct. It is less clear how one explicitly calculates a58 in algebraic form, except by brute force, and a considerable amount of brute force is required; but a numerical calculation to any reasonable accuracy is easily obtained from (7.3) and 1103 appears! The reader is encouraged to try this to, say, 16 digits. This presumably is what Ramanujan observed. Ironically, when Gosper computed 17 million digits of S using Sum 1, he had no mathematical proof that Sum 1 actually converged to 1/n. He compared ten million digits of the calculation to a previous calculation of Kanada et al. This verification that Sum 1 is correct to ten million places also provided the first complete proof that a58 is as advertised above. A nice touch-that the calculation of the sum should prove itself as it goes.
Roughly this works as follows. One knows enough about the exact algebraic nature of the components of dn(N) and XN to know that if the purported sum (of positive terms) were incorrect, that before one reached 3 million digits, this sum must have ceased to agree with 1/n. Notice that the components of Sum 1 are related to the solution of an equation of degree 58, but virtually no irrationality remains in the final packaging. Once again, there are very good number-theoretic reasons, presumably unknown to Ramanujan, why this must be so (58 is at least a good candidate number for such a reduction). Ramanujan's insight into this marvellous simplification remains obscure.
Ramanujan [26] gives 14 other series for 1/n, some others almost as spectacular as Sum 1-and one can indeed derive some even more spectacular related series.* He gives almost no explanation as to their genesis, saying only that there are "corresponding theories" to the standard theory (as sketched in section 5) from which they follow. Hardy, quoting Mordell, observed that "it is unfortunate that Ramanujan has not developed the corresponding theories." By methods analogous *(Added in proof) Many related series due to Borwein  This series is composed of fractions whose numerators grow like 26n and whose denominators are exactly 16-212n. In particular this can be used to calculate the second block of n binary digits of 7r without calculating the first n binary digits. This beautiful observation, due to Holloway, results, disappointingly, in no intrinsic reduction in complexity. 9. Sources. References [7], [11], [19], [26], [36], and [37] relate directly to Ramanujan's work. References [2], [8], [9], [10], [12], [22], [24], [27], [29], and [31] discuss the computational concerns of the paper.
Finally, details of all derivations are given in [11].