Hexadecimal Expansions

Last revised: September 9, 2003
Requires
Derivatives and Integrals
Power Series
References:
Pi Source Book, pages 663-676
Stan Wagon, Mathematica in Action, Chapter 20 (Science Center Library)
Mathematica proof of the BBP formula (attached)
Minilecture 1: Rewriting the arc tangent integral
Starting from

transform the integral to show that

The proof is elementary, although I have not seen it in print.  First make the substitution z = 1-y, which interchanges the limits of integration.  Then note that

All you have to do is multiply numerator and denominator by the appropriate factors to get
 in the denominator.
Minilecture 2: the original BBP integral
Starting from

transform the integral to show that

A derivation of this is in the attached Mathematica notebook file.
This integral can be done by partial fractions, though some help from Mathematica is nice.

To develop the original Bailey-Borwein-Plouffe integral

all you need to do is multiply numerator and denominator by 
This converts the integral to 

which is equivalent to equation (1.4) on p. 665 of the source book

Minilecture 3: The series expansion
Starting from 
p = , develop the Bailey-Borwein-Plouffe series

The secret is to expand
 = 
in a geometric series, then to integrate term by term.  Since
, you will get the crucial powers of 16 in the denominator.
Alternatively, from
p = 

you could develop the series expansion
,
which is mentioned on p. 665 of the Pi Source Book.

Minilecture 4: Extracting hex digits from the series
Explain how either of the above series can be used to compute a few consecutive hexadecimal digits of p, starting at position d in the hexadecimal expansion, using just ordinary computer double-precision floating-point arithmetic.
Once you have the series, using it to calculate the nth digit is still tricky.  The algorithm is explained on page 669 of the Pi Source Book and on pages 428 and 429 on Wagon's book, but it takes a while to understand why it doesn't require high-precision arithmetic.  You might try using the following steps to explain it.

1. To get the first three digits just after the decimal point in p, take
and multiply by 1000.  To get the sixth through eighth digits,
first multiply by : i.e. calculate

Note that 

2. So to get hex digits of p starting with the d+1st digit, we have to multiply p by and
take the result modulo 1.
With the series above this is easy because each term has a power of 16.
Take just the first term in the B-B-P series and multiply by , which leads to
 
The series now splits into two pieces, depending on whether or not there is a power of 16 in the denominator:

The second series (the infinite one) is easy, since it converges fast and the sum is less than 1.  If you want 5 hex digits after the decimal points, you won't need many more than 5 terms.
The first series looks hard, because it appears to involve large positive powers of 16.
3. But we only need to get the fractional part correct!
Consider the following:

You get the idea: we can take the numerator mod 7, and the fractional part is still the same.

4. Now there is an algorithm, attributed to Euclid, for efficiently raising an integer to the nth power.  First forget that the arithmetic is mod n.  The basic idea is illustrated by examples like
 16 in binary is (10000)
  31 in binary is (11111)
or  42 in binary is (101010)
The only operations required are squaring and multiplication by x.  Whether you multiply by x or not depends on whether a specfic bit in the binary expansion of the exponent is a 1 or a zero.  This algorithm is well specified in both references. 
5. But if we are going to divide an integer by (8i + 1) and keep just the fractional part, we can compute the integer using arithmetic mod (8i + 1) at every stage.  That means that we will never be squaring a number larger than 8d. Even if d is large, say a million, the square will be given correctly by double-precision floating point arithmetic
Minilecture 5: (computational):
Write a program that uses one or the other of the series above to calculate the 1000th through 1004th hexadecimal digits of p (where the 0th digit is 3, and the next 4 are 243F, so that p=3.243F...).  The correct answer, according to Stan Wagon(reference 2), is ...349F1.... 
Since this task doesn't require infinite-precision arithmetic, you can use C++ or
Mathematica.  (Mathematica has infinite-precision arithmetic, but you don't need it except to check your answer.)  Mathematica does have a built-in PowerMod function that does the most difficult part of the calculation for you.