Calculating A Billion Digits Of PI

One day I woke up and decided: I want to know as much PI as possible. And I mean "know" in the loosest way possible. As in, writing it down in a file.

I'm sure this has occurred to many people throughout human history. Lucky for me, I live in the modern age where I can read about anything I want, and I have metal and glass that can run computations for me. Sure enough, two days later I had a 1 GiB text file full of PI sitting on a VM I got from Azure. I could've gone much further, but I decided to stop then and there. It was a very fun two days, at some point the project stopped being about PI and became a problem of how to multiply really large integers efficiently. Let's talk about how I got there.

Bad

I started with what might be the silliest way to calculate PI, the sum from the Basel problem. This is a well known fact: $$\zeta(2) = \sum_{n=1}^{\infty} \dfrac1{n^2} = \dfrac{\pi^2}6$$

Rearrange and you've got PI. My code looked a little something like:

function pi_basel(nterms::Integer; bits::Integer = 256)
    setprecision(BigFloat, bits) do
        s = sum(inv(BigFloat(n)^2) for n in 1:nterms)
        return sqrt(6 * s)
    end
end

println(pi_basel(1_000))

After a good couple dozen minutes, I was maxxing out at only a few thousand digits. I was happy back then, but needless to say, this is a wonderfully useless solution. A similarly useless solution was the integral formulation: $$\int_0^1 \sqrt{1-x^2} \, dx = \dfrac\pi4$$

Next, I remembered Ramanujan has done some work on PI.

Not bad

Looking it up on Google, I came across this very nice note page. As the story goes, Ramanujan dreamt this up one night in typical Ramanujan fashion. Of course this was well grounded in his prior work and analysis. But I digress. The good news for me is that this formula is great: $$\dfrac1\pi = \dfrac{\sqrt{8}}{9801} \sum_{n=0}^\infty \dfrac{(4n)!}{(n!)^4} ​\dfrac{26390n + 1103}{396^{4n}}$$

It might as well be black magic to me! I know that this series inspired mathematicians to study a generalized form known as Ramanujan-Sato series, which are hypergeometric-like series resembling this one. But it produces a number of digits proportional to the number of terms you use, a linear relation. The number of digits per term works out to be around 8, genuinely great!

function pi_ramanujan(nterms::Integer; bits::Integer = 256)
    setprecision(BigFloat, bits) do
        factor = BigFloat(1) # Accumulated fraction factor
        total = BigFloat(0) # Accumulated inverse sum

        for n in 0:(nterms - 1)
            total += factor * (26390 * n + 1103)

            # Factor to accumulate (4n)!
            numerator =
                BigInt(4n + 1) *
                BigInt(4n + 2) *
                BigInt(4n + 3) *
                BigInt(4n + 4)

            # Factor to accumulate (n!)^4 and 396^4n
            denominator =
                BigInt(n + 1)^4 * BigInt(396)^4

            # Total factor to accumulate the fraction
            factor *= BigFloat(numerator) / BigFloat(denominator)
        end

        # Reciprocate and multiply by outside factors
        return BigFloat(9801) / (sqrt(BigFloat(8)) * total)
    end
end

println(pi_ramanujan(1_000))

A simple implementation like this performs remarkably well! In one go, I went from producing a few thousand digits straight to millions, on the same compute hardware (which would be my laptop).

The natural next thing to try is the Chudnovsky brothers formula, nicely noted in the page I linked above: $$\dfrac1\pi = \dfrac1{53360\sqrt{640320}} \sum_{n=0}^\infty (-1)^n \dfrac{(6n)!}{(n!)^3(3n)!} ​\dfrac{545140134n + 13591409}{640320^{3n}}$$

This series is due to the Chudnovsky brothers, and is the one commonly used to break PI digit world records. It asymptotically produces ~14 digits per term compared, significantly outperforming Ramanujan's formula.

I wrote some code for this one similar to the one above, and let it rip for about half an hour until I got 50 million digits. Not bad for a crude Julia implementation on a laptop!

Good

Throughout my research, I kept encountering the term "binary splitting". I found it on Wikipedia, random blog posts, published papers, and even Claude kept mentioning it when I asked it to review my code. Upon further inspection, I found out that binary splitting refers to a class of algorithms that efficiently evaluate series that can be shaped as a sum of products.

Before we get into what binary splitting is, we need to massage our series to take the form that binary splitting wants: $$\sum_n a_n \prod_k \dfrac{p_k}{q_k}$$

This involves quite a bit of algebraic gymnastics, so to keep this read compact, I've separated the proverbial elbow grease into its own page.

After that is done, we have our series in the binary splitting form: $$\sum_n a_n \prod_k \dfrac{p_k}{q_k} = \sum_n a_n \prod_k r_k$$

with: $$ a_n= 13591409+545140134n, $$

$$ p_k= -(6k-5)(2k-1)(6k-1), $$

$$ q_k= k^3\frac{640320^3}{24}. $$

Now let's talk about the binary splitting itself. Observe how each term of the sum has the same product as the previous sum term, but with one added product term.

Each term of the series as a column: a_n on the baseline with the factors r_1 through r_n stacked above it; every row is a single reused value∏k12345678∑n12345678a1r1a2r1r2a3r1r2r3a4r1r2r3r4a5r1r2r3r4r5a6r1r2r3r4r5r6a7r1r2r3r4r5r6r7a8r1r2r3r4r5r6r7r8×8×7×6×5×4×3×2×1all distinct

Each cell in this figure represents a product term. Multiplying the cells across one column, you get the corresponding sum term for that n value, with each row having a constant product index $k$.

An efficient implementation makes use of this by keeping track of the running product, and only multiplying it by the new product term for each new sum term. But binary splitting goes one step further. Start by separate the interval you're calculating on into two subintervals.

The range split in two: every column of the right block repeats the same factors r_1 to r_4, the product R(0,4) that the left block already forms∏k12345678∑n12345678a1r1a2r1r2a3r1r2r3a4r1r2r3r4a5r1r2r3r4r5a6r1r2r3r4r5r6a7r1r2r3r4r5r6r7a8r1r2r3r4r5r6r7r8R(0,4)S(0,4)S(4,8)

See the highlighted block in the second subinterval? Let's call it $X$. $X$ represents a constant product, repeated in every sum term, which happens to be the final value of the running product of the first subinterval. Here's the clever bit: Because it's repeated in every sum term, pull it out as a factor. You don't need to know what that value is to do the critical work of calculating the products above. This allows us to calculate the two subintervals in parallel, then merge them by multiplying the second subinterval's result by the left subinterval's $X$, then summing the two together.

Merging: the repeated factors collapse into the single product R(0,4), which the left block passes to the right block∏k12345678∑n12345678a1r1a2r1r2a3r1r2r3a4r1r2r3r4a5r5a6r5r6a7r5r6r7a8r5r6r7r8R(0,4)R(0,4)S(0,4), R(0,4)S(4,8), R(4,8)

That was a mouthful, so let's formulate it into notation. We'll look at the general case of a subinterval $[u, v)$. We want our subinterval calculation to result with two characteristic values:

  • The running product, which we called $X$ prior. $$R(u, v) = \prod_{k=u+1}^{v} r_k = \prod_{k=u+1}^{v} \dfrac{p_k}{q_k}$$
  • The running sum, which doesn't include the product terms we decided to factor out. $$S(u, v) = \sum_{n=u+1}^{v} a_n \prod_{k=u+1}^{n} r_k$$

Now, of course, each subinterval will internally calculate these two incrementally, with the running product receiving an additional term each time. Now, we can write our merging rule in notation. Given $a < b < c$, we deduce: $$S(a, c) = S(a, b) + R(a, b) S(b, c)$$

This pops out from the formulae above, but it is exactly the intuition shown in the diagram above, factoring out the block that's already calculated by the first subinterval.

We have effectively separated the calculation of subintervals, and achieved our binary splitting. But why stop at one splitting? You can split your desired full interval as many times as you'd like, untill the intervals are small enough to trivially calculate. The great part is that you can calculate the subintervals and merge them in parallel, making efficient use of your CPU cores. This is precisely how world records are achieved.

Best

Well, not really the best, but the best I did! There is one last improvement we can make. Notice how the running product $R$ is a fraction, so implementing this with BigRational datatypes results in a huge performance cost every time a merge runs. The addition becomes very costly, having to run rational reduction and denominator alignment on extremely large numerator and denominator fractions. What we can do is massage this formulation a little further to make it pure integers only.

We define two auxiliary integers, which are essentially the numerator and denominator versions of $R$: $$P(u, v) = \prod_{k=u+1}^{v} p_k$$ $$Q(u, v) = \prod_{k=u+1}^{v} q_k$$

Then notice how in its unreduced form, $S(u, v)$ has a denominator of $Q(u, v)$. Well then, if we define our final auxiliary function: $$T(u, v) = S(u, v) Q(u, v)$$

We know for sure it is guaranteed to be an integer. Now to derive the splitting formula for it. We know: $$S(a, c) = S(a, b) + \dfrac{P(a, b)}{Q(a, b)} S(b, c)$$

Multiply by $Q(a, c) = Q(a, b) Q(b, c)$: $$Q(a, c) S(a, c) = Q(b, c) Q(a, b) S(a, b) + P(a, b) Q(b, c) S(b, c)$$

Sub in the $T$ definition: $$T(a, c) = Q(b, c) T(a, b) + P(a, b) T(b, c)$$

So now, our complete merge rule is: $$Q(a, c) = Q(a, b) Q(b, c)$$ $$P(a, c) = P(a, b) P(b, c)$$ $$T(a, c) = Q(b, c) T(a, b) + P(a, b) T(b, c)$$

And everything is an integer!

Some further reduction optimizations are possible if you're clever about it, but this is where I stopped. I highly recommend the technical paper by Bellard for a much more formal analysis and for those further optimizations.

Here's a clean implementation in Julia that does single-core splitting:

const CHUD_A = BigInt(13_591_409)
const CHUD_B = BigInt(545_140_134)
const CHUD_C = BigInt(640_320)

const C3_OVER_24 = CHUD_C^3 ÷ 24

function binary_split(left::Int, right::Int)::NTuple{3, BigInt}
    # Leaf of the recursion tree.
    if right - left == 1
        if left == 0
            return (BigInt(1), BigInt(1), CHUD_A)
        end

        n = BigInt(left)

        P = (6*n - 5) * (2*n - 1) * (6*n - 1)
        Q = -n^3 * C3_OVER_24
        T = P * (CHUD_A + CHUD_B*n)

        return (P, Q, T)
    end

    # Split the interval in half.
    middle = (left + right) ÷ 2

    Pl, Ql, Tl = binary_split(left, middle)
    Pr, Qr, Tr = binary_split(middle, right)

    # Merge the two adjacent intervals.
    P = Pl * Pr
    Q = Ql * Qr
    T = Tl * Qr + Pl * Tr

    return (P, Q, T)
end

Turning the result into PI is a matter of calling with the desired number of terms, setting the appropriate floating point precision, folding in the constants, and converting to decimal.

Ugly

Here's my incredibly cursed full implementation that I used to achieve a billion digits:

println("Starting...")

target_digits = 1_000_000_000
N = BigInt(ceil(target_digits / 14.18))
bits = Int(ceil(target_digits * 1.1 * log2(10)))
setprecision(BigFloat, bits)

println("Calculating PI to $target_digits digits, using $N terms of the Chudnovsky formula and $bits bits of precision...")

const A = BigInt(545140134)
const B = BigInt(13591409)
const Z = (BigInt(640320)^3) ÷ 24

function calculate_split(left::BigInt, right::BigInt)::Tuple{BigInt, BigInt, BigInt}
    if right - left == 1
        if left == 0
            return (BigInt(1), BigInt(1), BigInt(B))
        else
            P = (6 * left - 5) * (6 * left - 1) * (2 * left - 1)
            Q = -left^3 * Z
            T = P * (A * left + B)
            return (P, Q, T)
        end
    else
        middle = (left + right) ÷ 2
        Pl, Ql, Tl = calculate_split(left, middle)
        Pr, Qr, Tr = calculate_split(middle, right)

        P = Pl * Pr
        Q = Ql * Qr
        T = Tl * Qr + Tr * Pl

        return (P, Q, T)
    end
end

if Threads.nthreads() < 8
    println("The number of threads is less than 8. Cannot cleanly split the work.")
    exit();
end

interval_size = BigInt(N) ÷ BigInt(16)
PQTs = Array{Tuple{BigInt, BigInt, BigInt}}(undef, 16)

println("Initial calculation starting...")

function run_thread(i::Int)
    left_left = i * interval_size
    left_right = (i + 1) * interval_size
    right_left = (15 - i) * interval_size
    right_right = (16 - i) * interval_size

    Pl, Ql, Tl = calculate_split(left_left, left_right)
    Pr, Qr, Tr = calculate_split(right_left, right_right)
    PQTs[i + 1] = (Pl, Ql, Tl)
    PQTs[16 - i] = (Pr, Qr, Tr)
end

@time Threads.@threads for i in 0:7
    run_thread(i)
end

println("Merging the results from the threads...")

# Recursively merge the results from the threads
function merge_results(PQTs::Array{Tuple{BigInt, BigInt, BigInt}})
    # Construct an array half the size of PQTs to perform one merge step
    # Recurse with this function until we have the top of the tree
    if length(PQTs) == 4
        merged = Array{Tuple{BigInt, BigInt, BigInt}}(undef, length(PQTs) ÷ 2)
        Threads.@threads for i in 1:length(merged)
            Pl, Ql, Tl = PQTs[2 * i - 1]
            Pr, Qr, Tr = PQTs[2 * i]

            P = BigInt(0)
            Q = BigInt(0)

            P = 0
            if i == 1
                P = Pl * Pr
            end
            Q = Ql * Qr

            T = Tl * Qr + Tr * Pl
            merged[i] = (P, Q, T)
        end

        return merged
    else
        merged = Array{Tuple{BigInt, BigInt, BigInt}}(undef, length(PQTs) ÷ 2)
        Threads.@threads for i in 1:length(merged)
            Pl, Ql, Tl = PQTs[2 * i - 1]
            Pr, Qr, Tr = PQTs[2 * i]
            P = Pl * Pr
            Q = Ql * Qr
            T = Tl * Qr + Tr * Pl
            merged[i] = (P, Q, T)
        end

        return merge_results(merged)
    end
end

@time begin
    (Pl, Ql, Tl), (_, Qr, Tr) = merge_results(PQTs)
    final_results = Array{BigInt}(undef, 3)

    Threads.@threads for i in 1:3
        if i == 1
            final_results[1] = Tl * Qr
        elseif i == 2
            final_results[2] = Tr * Pl
        else
            final_results[3] = Ql * Qr
        end
    end

    T = final_results[1] + final_results[2]
end

println("Calculating the floating-point result...")
c = BigFloat(53360) * sqrt(BigFloat(640320))
@time calculated_pi = c * BigFloat(final_results[3]) / BigFloat(T)

println("Saving raw result to txt file...")
open("pi.txt", "w") do f
    @time println(f, calculated_pi)
end