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^-n 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^n-1, 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 ofx. 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 change 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 calculations 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 we did not need 100 hex digits (400 bits) of precision to calculate the 100th hex digit of . We used standard precision, which is between 15 and 16 bits.

Improved code

We can improve the code above by adding an estimate of the infinite series we ignored.

A more subtle improvement is reducing the sum accumulator variable s mod 1. We only need the fractional part of s, and so by routinely cutting off the integer part we keep s from getting large. This improves accuracy by devoting all the bits in the machine representation of s to the fractional part.

epsilon = np.finfo(float).epsdef term(N, n): return 16**(N-1-n) * (4/(8*n + 1) - 2/(8*n+4) - 1/(8*n+5) - 1/(8*n + 6))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) s = s%1 n = N + 1 while True: t = term(N, n) s += t n += 1 if abs(t) < epsilon: break frac = s - floor(s) return floor(16*frac)

I checked the output of this code with [1] for N = 10⁶, 10⁷, 10⁸, and 10⁹. The code will fail due to rounding error for some very large N. We could refine the code to push this value of N a little further out, but eventually machine precision will be the limiting factor.

[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.

The post Computing the nth digit of directly first appeared on John D. Cook.