Computing the nth digit of π directly

The following equation, known as the BBP formula [1], will let you compute the nth digit of π directly without having to compute the previous digits.

I’ve seen this claim many times, but I’ve never seen it explained in much detail how this formula lets you extract digits of π.

First of all, this formula lets you calculate the digits of π, not in base 10, but in base 16, i.e. in hexadecimal.

It looks like the BBP formula might let you extract hexadecimal digits. After all, the hexadecimal expansion of π is the set of coefficients a_n such that

where each a_n is an integer between 0 and 15. But the term multiplying 16⁻ⁿ in the BBP formula is not an integer, so it’s not that easy. Fortunately, it’s not that hard either.

Warmup

As a warmup, let’s think how we would find the nth digit of π in base 10. We’d multiply by 10ⁿ, throw away the factional part, and look at the last digit. That is, we would calculate

Another approach would be to multiply by 10ⁿ⁻¹, shifting the digit that we’re after to the first digit after the decimal point, then take the integer part of 10 times the fraction part.

Here {x} denotes the fractional part of x. This is a little more complicated, but now all the calculation is inside the floor operator.

For example, suppose we wanted to find the 4th digit of π. Multiplying by 10³ gives 3141.592… with fractional part 0.592…. Multiplying by 10 gives 5.92… and taking the floor gives us 5.

100th hex digit

Now to find the nth hexadecimal digit we’d do the analogous procedure, replacing 10 with 16. To make things concrete, let’s calculate the 100th hexadecimal digit. We need to calculate

We can replace the infinite sum with a sum up to 99 because the remaining terms sum to an amount too small to effect our answer. Note that we’re being sloppy here, but this step is justified in this particular example.

Here’s the trick that makes this computation practical: when calculating the fractional part, we can carry out the calculation of the first term mod (8n + 1), and the second part mod (8n + 4), etc. We can use fast modular exponentiation, the same trick that makes, for example, a lot of encryption calculatons practical.

Here’s code that evaluates the Nth hexadecimal digit of π by evaluating the expression above.

def hexdigit(N):
    s = 0
    for n in range(N):
        s += 4*pow(16, N-n-1, 8*n + 1) / (8*n + 1)
        s -= 2*pow(16, N-n-1, 8*n + 4) / (8*n + 4)
        s -=   pow(16, N-n-1, 8*n + 5) / (8*n + 5)
        s -=   pow(16, N-n-1, 8*n + 6) / (8*n + 6)
    frac = s - floor(s)
    return floor(16*frac)

Here the three-argument version of pow, introduced into Python a few years ago, carries out modular exponentiation efficiently. That is, pow(b, x, m) calculates b^x mod m.

This code correctly calculates that the 100th hex digit of π is C. Note that because we did not include an estimate of the tail of the sum, this code could fail for some values of N. A more rigorous version of the code would add a tail estimate to the variable s above.

[1] Named after the initials of the authors. Bailey, David H.; Borwein, Peter B.; Plouffe, Simon (1997). On the Rapid Computation of Various Polylogarithmic Constants. Mathematics of Computation. 66 (218) 903–913.