Read the latest issue of Notices  Read the latest issue of Bulletin  Shop in the AMS Bookstore  My Account | Cart  
 
American Mathematical Society   

Pulling Digits out of Pi

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

Email to a friendMail to a friend Print this articlePrint this article

 

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 50-digit 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 an denotes the circumference of the circumscribed n-gon (shown in red) and bn the circumference of the inscribed n-gon (shown in blue), we have

In the case that n = 6, the circumferences are easily found to be

Essentially using half-angle 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 = -u2, 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 high-precision 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 high-precision 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 dn+1, we multiply by 2n:

Notice that the digit we wish to find, dn+1, is the first digit behind the decimal point. Therefore, dn+1 = 0 if the fractional part of 2n ln(2) is smaller than 0.5. Otherwise, dn+1 = 1. Therefore, we simply need to compute the fractional part of 2n ln(2), denoted { 2n ln(2) }, which we do as follows:

The terms in the first sum may be computed efficiently using the well-known 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 16k 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 (1015) 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 year-long 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 so-called 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 x1, x2, ... xn, we may ask if there are integers a1, a2, ... an, at least one of which is not zero, such that

a1x1 + a2x2 + ... anxn = 0.

For a given set of real numbers, there may be no set of integers a1, a2, ... an 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 a1 and a2 such that a1 A + a2 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 well-known 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 a1 and a2.

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 round-off 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 nth 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

Integer Relation Detection

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.