# Computing the hexadecimal value of pi

The
Bailey–Borwein–Plouffe formula is
one of the
several
algorithms to compute $\pi$.
Here it is:

What makes this formula stand out among other approximations of $\pi$ is that
it allows one to directly extract the $n$-th fractional digit of the
hexadecimal value of $\pi$ without computing the preceding ones.

Image credit: Cormullion, Julia
code
here.

The Wikipedia article about the Bailey–Borwein–Plouffe formula explains that the
$n$-th fractional digit (well, actually it’s the $n+1$-th) $d_{n}$ is
given by

where

Only the fractional part of expression in square brackets on the right side of
$d_{n}$ is relevant, thus, in order to avoid rounding errors, when we compute
each term of the finite sum above we can take only the fractional part. This
allows us to always use ordinary double precision floating-point arithmetic,
without resorting to arbitrary-precision numbers. In addition note that the
terms of the infinite sum get quickly very small, so we can stop the summation
when they become negligible.

## Implementation in Julia language

Here is a Julia implementation of the algorithm to
extract the $n$-th fractional digit of $\pi$:

The pi_digit function gives the $n$-th hexadecimal fractional digit of
$\pi$ as a base-10 integer, and the pi_string function returns the first
$n$ hexadecimal digits of $\pi$ as a valid hexadecimal floating-point
literal:

julia> pi_digit(1)
2

julia> pi_digit(6)
10

julia> pi_string(1000)


While I was preparing this post I found an unregistered
package PiBBP.jl that implements the
Bailey–Borwein–Plouffe formula. This is faster than my code above, mostly
thanks to a function
for
modular exponentiation
more efficient than that available in Julia standard library.

## Test results

Let’s check if the function is working correctly. We can use
the
parse function
to convert the string to a decimal floating point
number. IEEE 754 double precision
floating-point numbers have a 53-bit mantissa, amounting to $53 / \log_{2}(16) \approx 13$ hexadecimal digits:

julia> pi_string(13)
"0x3.243f6a8885a30p0"

julia> parse(Float64, pi_string(13))
3.141592653589793

julia> Float64(π) == parse(Float64, pi_string(13))
true


Generator expressions allow
us to obtain the decimal value of the number in a very simple way, without using

julia> 3 + sum(pi_digit(n)/16^n for n in 1:13)
3.141592653589793


We can use
the
arbitrary-precision BigFloat
to check the correctness of the result for even more digits. By default,
BigFloat numbers in Julia have a 256-bit mantissa:

julia> precision(BigFloat)
256


The result is correct for the first $256 / \log_{2}(16) = 64$ hexadecimal
digits:

julia> pi_string(64)
"0x3.243f6a8885a308d313198a2e03707344a4093822299f31d0082efa98ec4e6c89p0"

julia> BigFloat(π) == parse(BigFloat, pi_string(64))
true

julia> 3 + sum(pi_digit(n)/big(16)^n for n in 1:64)
3.141592653589793238462643383279502884197169399375105820974944592307816406286198


It’s possible to increase the precision of BigFloat numbers, to further test
the accuracy of the Bailey–Borwein–Plouffe formula:

julia> setprecision(BigFloat, 4000) do
BigFloat(π) == parse(BigFloat, pi_string(1000))
end
true


Since the Bailey–Borwein–Plouffe formula extracts the $n$-th digit of $\pi$
without computing the other ones, we can write a multi-threaded version of
pi_string, taking advantage of native support
for
in Julia:

For example, running Julia with 4 threads gives a 2× speed-up:

julia> Threads.nthreads()
4

julia> using BenchmarkTools

true

julia> @benchmark pi_string(1000)
BenchmarkTools.Trial:
memory estimate:  105.28 KiB
allocs estimate:  2016
--------------
minimum time:     556.228 ms (0.00% GC)
median time:      559.198 ms (0.00% GC)
mean time:        559.579 ms (0.00% GC)
maximum time:     564.502 ms (0.00% GC)
--------------
samples:          9
evals/sample:     1