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.