Pulling Digits out of Pi
Posted December 2007.
In this article, we will investigate how , the ratio of the circumference of a circle to its diameter, may be computed to an astonishing degree of accuracy...
David Austin
Grand Valley State University
david at merganser.math.gvsu.edu
Introduction
Nature gives us many important numbers, such as and e, that are not easily represented by human number systems. For instance, the first fifty digits of
= 3.14159265358979323846264338327950288419716939937510...
provide only an approximation of this important number, and we feel compelled to understand numbers like this more deeply. We might also ask if there is any relationship between important numbers. For instance, at first glance, and e seem to arise from different areas of mathematics so it would be remarkable to discover that their numerical values are related in some simple way.
In this article, we will investigate how , the ratio of the circumference of a circle to its diameter, may be computed to an astonishing degree of accuracy. In fact, we'll see how to reach deeply inside the hexadecimal representation of to compute any digit that we wish. We'll also describe a related new algorithm that allows us to investigate whether simple relationships exist between various numerical constants.
Of course, we would never need to know with 50digit accuracy for any "real world" application. The history of computing , however, is rich so these ideas continue a venerable current of mathematical investigation motivated mainly by our own curiosity. For computational scientists, computing presents technical challenges, such as how to multiply numbers efficiently, whose solutions are broadly useful.
Earlier computations of pi
The Greek mathematician Archimedes was one of the first to undertake a careful study of . His technique begins with the unit circle, whose circumference is , and approximates its circumference by the perimeters of inscribed and circumscribed regular polygons.
For instance, the perimeters of the two hexagons shown below will lead to lower and upper bounds on the value of .
We obtain better bounds by increasing the number of sides of the polygons.
If a_{n} denotes the circumference of the circumscribed ngon (shown in red) and b_{n} the circumference of the inscribed ngon (shown in blue), we have
In the case that n = 6, the circumferences are easily found to be
Essentially using halfangle identities from trigonometry, Archimedes determined how the perimeter changes when the number of sides doubles to obtain the recurrence relations
In this way, he found the bounds
Ludolph van Ceulen, a German mathematician working in the Netherlands around 1600, used Archimedes' technique to compute the first 35 digits of and had the lower and upper approximations inscribed on his tombstone.

Photograph courtesy of Karen Aardal 
Not long after, the development of calculus gave new means of computing , one of which we'll see now. Remember that a geometric series has the form:
That is, each term in the sum is formed by multiplying its predecessor by a fixed number r. If the absolute value of r is less than one, the series converges and may be easily evaluated. To see this, call the sum S and multiply by r as below:
Now substract the second row from the first to obtain (1  r) S = 1 so that
If we let r = u^{2}, we obtain
and then
If we set x = 1, we then obtain an expression for that is often attributed to Leibniz:
In practice, this is not a particularly useful way to compute since it converges so slowly. Here are the approximations given by the first few terms:
n 
Sum of first n terms 
1 
4.0000 
2 
2.6667 
3 
3.4667 
4 
2.8952 
5 
3.3397 
6 
2.9760 
... 
... 
1000 
3.1406 
1001 
3.1426 
As you can see, after adding a thousand terms, we have only found three digits of .
Euler found a more useful formula,
, 
which may be simply explained using complex multiplication. Remember that the argument of the complex product ab is the sum of the arguments of a and b. Euler's formula then results from the fact that (2 + i) (3 + i) = 5 + 5i. Using our series for the arctangent function, we are led to
Adding together the first ten terms of both series gives the approximation 3.1415925796, accurate in the first seven digits. This, and other similar formulas involving the arctangent, eventually produced the first few hundreds of digits of .
Of course, all of this was accomplished before the advent of mechanical calculation machines. Recently, Yasumasa Kanada has computed over one trillion digits of using a similar arctangent formula and about 600 hours of processing time on a supercomputer.
The BBP formula
If we would like to compute a particular digit of , the techniques outlined above require us to compute all the digits that come before using highprecision arithmetic. In the mid 1990's, a remarkable new formula for was discovered by David Bailey, Peter Borwein, and Simon Plouffe (BBP):
We'll describe the discovery of this formula a bit later. First, we'll see how this allows us to compute a hexadecimal digit directly without computing all the digits that come before it and without using highprecision arithmetic.
To understand how this works, let's consider a similar type of formula for ln(2)=0.69314718056.... Remember that
so that
Let's now look at the binary digits of ln(2):
If we want to compute the (n+1)^{th} digit d_{n+1}, we multiply by 2^{n}:
Notice that the digit we wish to find, d_{n+1}, is the first digit behind the decimal point. Therefore, d_{n+1} = 0 if the fractional part of 2^{n} ln(2) is smaller than 0.5. Otherwise, d_{n+1} = 1. Therefore, we simply need to compute the fractional part of 2^{n} ln(2), denoted { 2^{n} ln(2) }, which we do as follows:
The terms in the first sum may be computed efficiently using the wellknown binary exponentiation algorithm. The second term is smaller than 1/(n+1) so if n is large, the second term makes only a small contribution. Remembering that we only need to determine whether this fractional part is smaller or larger than 0.5, we simply compute enough terms in the second sum to settle the question.
Turning back to the BBP formula
we see that the presence of 16^{k} in each of the terms allows us to compute, in a similar fashion, chosen hexadecimal digits of without computing the intermediate digits.
The hexadecimal digits, beginning at the trillion and first, are in fact B4466E8D215388C4E014. Recent work has succeeded in computing the quadrillionth (10^{15}) hexadecimal digit of .
The PSLQ Algorithm
One way to prove the BBP formula is by considering the integral
On one hand, the technique of partial fractions, familiar to most students who have taken a yearlong calculus course, shows that this integral equals . On the other, the geometric series shows that
a fact that enables us to write the integrand as a series that may be integrated term by term. Putting these two facts together gives the BBP formula.
But how was the BBP formula discovered? Certainly not through the integral just described. Rather, the BBP formula results from the socalled PSLQ algorithm, an algorithm designed to detect integer relations between sets of real numbers.
For instance, if we are given a set of real numbers x_{1}, x_{2}, ... x_{n}, we may ask if there are integers a_{1}, a_{2}, ... a_{n}, at least one of which is not zero, such that
a_{1}x_{1} + a_{2}x_{2} + ... a_{n}x_{n} = 0.
For a given set of real numbers, there may be no set of integers a_{1}, a_{2}, ... a_{n} giving a relation as above. If there is, however, the PSLQ algorithm will find one such set.
Detecting integer relations has a long history. For example, Euclid considered this problem in Book X of The Elements for the case when n = 2. To explain Euclid's algorithm, we will call the real numbers A and B. The requirement that there be integers a_{1} and a_{2} such that a_{1} A + a_{2} B = 0 is the same as asking if there is a real number C and integers n and m such that A = nC and B = mC. (In this case, Euclid said that A and B are commensurable.)
Euclid's wellknown algorithm works in the following way for positive A and B.
 If A = B, then there is an integer relation between A and B and the algorithm terminates with C = A = B.
 Call B the smaller of the two real numbers and let A = A  B.
 Repeat.
If the algorithm terminates, then we have found C so that A = nC and B = mC. If we let
thus demonstrating an integer relation.
Two illustrations, in which the lengths of the bars represent A and B, of the algorithm are shown below. In the first case, we find that A and B are commensurable and the value of C is indicated.
In the next case, A = and B = e.
After several steps, the algorithm has not terminated. Of course, this does not mean that it won't if we let it run a while longer; in practice, we cannot determine that the algorithm will run indefinitely. We can, however, note that if the algorithm has not terminated, then we can guarantee that C is smaller than the current value of A. This in turn gives us lower bounds on the integers a_{1} and a_{2}.
Since Euclid's time, the problem of generalizing Euclid's algorithm to sets of more than two real numbers was considered by Euler, Jacobi, and Poincare, among others. In 1979, an algorithm was found by Helaman Ferguson and Rodney Forcade and subsequently refined to the PSLQ algorithm by Ferguson, Bailey and Steve Arno.
The PSLQ algorithm is similar to the Euclidean algorithm, which it generalizes: If there is an integer relation, the algorithm will terminate when it is found. If the algorithm does not terminate, there is no integer relation. Of course, we can never know in practice that the algorithm does not terminate; however, if the algorithm has not terminated after running for a while, we can determine lower bounds on the size of the integers in a relation.
Implementing the PSLQ algorithm presents its own challenges since computers can only work with finite precision arithmetic. For instance, the real numbers for whom we seek an integer relation can only be approximated by finitely many digits in computer memory. In addition, the real arithmetic needed in the algorithm will be subject to roundoff error.
Because of this, implementations of the PSLQ algorithm are not able to produce definite integer relations. Instead, the algorithm suggests likely relations, the proofs of which need to be provided independent of the algorithm.
The BBP formula was discovered using the PSLQ algorithm, along with what the authors call "inspired guesswork and extensive searching," to find an integer relation between and the constants
where j = 1, 2, 3, ..., 8. The formula was then proven using the integral argument given above.
The BBP formula enables us to compute the n^{th} hexadecimal digit of without computing the first n  1 digits. At this time, it is not known if a similar formula exists that would allow us to compute arbitrary decimal digits in the same way.
Besides the discovery of the BBP formula for , the PSLQ algorithm has found other uses. For instance, the algorithm may be used to investigate whether a given real number is an algebraic number (that is, a root of a polynomial with integer coefficients). In addition, PSLQ has found identities involving constants that arise in quantum field theory from the evaluation of certain Feynman diagrams.
Incidentally, Helaman Ferguson, whose mathematical work has been fundamental in developing integer relation detecting algorithms, is also a sculptor whose art expresses both beauty and mathematical understanding.
References
Survey articles
 David H. Bailey and Jonathan M. Borwein, Experimental Mathematics: Examples, Methods and Implications, Notices of the American Mathematical Society, Vol. 52 (5), 2005, 502  514.
 Barry Cipra, Digits of Pi, What's Happening in the Mathematical Sciences, Vol. 6, 2006, 28  39.
Integer Relation Detection
 H.R.P. Ferguson, R.W. Forcade, Generalization of the Euclidean Algorithm for real numbers to all dimensions higher than two, Bulletin of the American Mathematical Society, Vol. 1 (6), 1979, 912  914.
 Helaman R.P. Ferguson, David H. Bailey, Steve Arno, Analysis of PSLQ, an integer relation finding algorithm, Mathematics of Computation, Vol. 68 (225), 1999, 351369.
 David H. Bailey, Helaman R.P. Ferguson, Numerical results on relations between fundamental constants using a new algorithm, Mathematics of Computation, Vol. 53 (188), 1989, 649  656.
 David H. Bailey, David J. Broadhurst, Parallel integer relation detection: techniques and applications, Mathematics of Computation, Vol. 70 (236), 1719  1736.
The BPP paper
The Euclidean algorithm
David Austin
Grand Valley State University
david at merganser.math.gvsu.edu
Those who can access JSTOR can find some of the papers mentioned above there. For those with access, the American Mathematical Society's MathSciNet can be used to get additional bibliographic information and reviews of some these materials. Some of the items above can be accessed via the ACM Portal , which also provides bibliographic services.