Algorithm for Calculating Individual Hexadecimal Digits of Pi
The algorithm for computing individual hex digits of pi is based on
this formula:
infinity
----
\ 1 [ 4 2 1 1 ]
pi = | ---- [---- - ---- - ---- - ----]
/ 16^k [8k+1 8k+4 8k+5 8k+6]
----
k=0
The algorithm is as follows: Write the above formula as four
individual sums, and let s1 be the first sum, without the 4:
infinity
----
\ 1
s1 = | ----------
/ 16^k(8k+1)
----
k=0
The hex (base 16) expansion of s1 beginning at position d+1 is given
by frac(16^d s1), where frac() denotes fractional part. We can write
[ d ] [infinity ]
[---- ] [ ---- ]
[\ 16^(d-k)] [ \ 16^(d-k)]
frac(16^d s1) = frac[| --------] + frac[ | --------]
[/ 8k+1 ] [ / 8k+1 ]
[---- ] [---- ]
[ k=0 [k=d+1 ]
[ d ] [infinity ]
[---- ] [ ---- ]
[\ 16^(d-k) mod 8k+1] [ \ 16^(d-k)]
= frac[| -----------------] + frac[ | --------]
[/ 8k+1 ] [ / 8k+1 ]
[---- ] [---- ]
[ k=0 [k=d+1 ]
Note that we are justified in inserting "mod 8k+1" in the numerator of
the first sum, since we are only interested in the fractional part of
the result.
The key "trick" is this: the numerator of the first sum can be very
rapidly calculated by employing the binary algorithm for
exponentiation, where we reduce modulo 8k+1 after each multiplication.
As an illustration, we can compute 3^16 mod 100 as follows:
3^2 = 9; 9 mod 100 = 9
9^2 = 81; 81 mod 100 = 81
81^2 = 6561; 6561 mod 100 = 61
61^2 = 3721; 3721 mod 100 = 21
As a check, 3^16 = 43,046,721, so that 3^16 mod 100 is indeed 21.
Note that the above calculation was performed using only 4
muliplications and 4 mod operations, with integers that do not exceed
100^2 = 10,000 in size. In the pi algorithm we are describing,
ordinary 64-bit integer arithmetic suffices to calculate the
375 millionth hex digit of pi.
After each numerator is calculated using this scheme, the result is
divided by 8k+1 using ordinary IEEE floating-point arithmetic, and
then added to the sum, discarding integer parts whenever they appear.
For the second sum, only a few terms need be calculated, using IEEE
floating-point arithmetic, since it rapidly converges. The final
result, when expressed as a hexadecimal number, contains the desired
hex expansion of s1 beginning at position d+1. Doing this for each of
the four sums s1, s2, s3, s4 gives the hex expansion of pi beginning
at position d+1, plus about 8 or so additional hex digits.
C and Fortran-90 programs are available that implement this scheme. See
the main page for links to these files.