How to calculate 1 million digits of pi

\(\pi\) day (half \(\tau\) day? hmmm) rolled around last week, and it got me wondering - how exactly does one go about calculating large numbers of digits for mathematical constants, \(\pi\) in particular?  One hears about how \(\pi\) , \(e\) , and the like have been calculated out to millions, billions, or even trillions of digits.  I had a vague idea that some series expansion is used in order to build up the value piece by piece.  But is each digit calculated individually, or does the entire value need to be stored in memory at once?  What kind of data types are needed?

Some simple web searches on the subject will steer you toward the Wikipedia articles on \(\pi\)  and approximations of \(\pi\) .  These give nice examples of the various different series expansions which can be used, as well as some interesting history on the topic.  We note that the Baily-Borwein-Plouffe formula or its relatives can be used to calculate (base-16) digits of \(\pi\) independently (i.e. without needing to calculate any other digits besides the one you want), but that this is somewhat of a special case.  There is no general method for doing this for an arbitrary mathematical constant.  Until someone comes up with a specialized clever formula, the entire value up to the desired precision typically needs to be calculated together.

What’s missing from the Wikipedia articles are implementation details.  Further web searches took me to a few sites with code that was fabulously confusing, or even plain wrong.  However this site, last updated in 1997, turned out to be just what I was looking for.  The author provides some of the details of how to do the calculation via the quickly-converging formula (found by Machin in 1706)

He gives a great explanation of how to code the Taylor expansion of arctan, and  various pitfalls of different computational approaches.  However, even this page leaves out some crucial details.  He refers to the formula above, as well as the Taylor expansion itself, as “black magic” which should be taken on faith.  That’s no good if you ask me.  I don’t know how Machin came up with his formula, but it is not too hard to verify that it’s accurate.  And many of us were taught Taylor series in high school, so they can’t be all that bad. [Historical bonus facts: Machin was Taylor’s doctoral advisor at Cambridge, and the two served together in the Royal Society, both ending up on the board which adjudicated the Newton/Leibniz dispute.]

The biggest omission, though, is that the author doesn’t reveal how to implement the high-precision number class he’s using in all of his code.  He teases us by merely stating that with “careful coding of these routines” we can do the calculation.  He also makes frequent reference to C++ being a “high-level” language, which I guess was less amusing of a statement in 1997.

So let’s create our own high-precision number class (in C#) and carry out the calculation.  Here are the goals:

1. Use only types that are in standard .NET libraries (.NET 4.0 is most recent release at the time of writing)
   • OK I’m cheating a bit here.  We will use the unbounded integer class System.Numerics.BigInteger, which is fantastically optimized and will give us a lot for free.  Trust me, though, our code would be fast enough even with a super simple BigDecimal implementation (i.e. just keep the decimal digits before and after the decimal point in int arrays, and carry out normal long division, multiplication, etc, like we learned in elementary school).  This was actually my first approach, but it requires a lot more code.
2. Have it be fast enough to calculate 1 million digits in “reasonable” time.

Our HighPrecision class will be a limited-precision rational number class, which only keeps track of enough data to maintain some set number of accurate digits after the decimal.  How do we do that?  Well say we some rational number \(r = \frac{p}{q}\) and we want to find its decimal representation to \(k\) digits after the decimal.  Then what we really want is some rational number \(r^{\prime} = a/10^k\) , since any decimal number with a finite number of digits can be expressed as an integer divided by an appropriately large power of 10 (e.g. 0.25 = ²⁵⁄₁₀₀).   To calculate \(a\) , we simply set \(a/10^k=\frac{p}{q}\) and come to \(a = p \times 10^k / q\) .  We can let this quotient be truncated like normal with computer integer division, since we are only hoping for an approximation to \(r\) .

So really, the trick to our limited-precision rational class is to keep a constant denominator equal to \(10^k\) , and make sure that during any calculation the numerator is appropriately rescaled.  Here’s the code:

public class HighPrecision
{
    private static BigInteger denom;
    private static int precision;
    private static int slop = 4;
    private BigInteger num;

    public HighPrecision(BigInteger numerator, BigInteger denominator)
    {
        // public constructor rescales numerator as needed
        num = (HighPrecision.denom * numerator) / denominator;
    }

    private HighPrecision(BigInteger numerator)
    {
        // private constructor for when we already know the scaling
        num = numerator;
    }

    public static int Precision
    {
        get { return precision; }
        set
        {
            HighPrecision.precision = value;
            denom = BigInteger.Pow(10, HighPrecision.precision + slop);  // leave some buffer
        }
    }

    public bool IsZero
    {
        get { return num.IsZero; }
    }

    public BigInteger Numerator
    {
        get { return num; }
    }

    public BigInteger Denominator
    {
        get { return HighPrecision.denom; }
    }

    public static HighPrecision operator *(int left, HighPrecision right)
    {
        // no need to resale when multiplying by an int
        return new HighPrecision(right.num * left);
    }

    public static HighPrecision operator *(HighPrecision left, HighPrecision right)
    {
        // a/D * b/D = ab/D^2 = (ab/D)/D
        return new HighPrecision((left.num * right.num) / HighPrecision.denom);
    }

    public static HighPrecision operator /(HighPrecision left, int right)
    {
        // no need to rescale when dividing by an int
        return new HighPrecision(left.num / right);
    }

    public static HighPrecision operator +(HighPrecision left, HighPrecision right)
    {
        // when we know the denominators are the same, can just add the numerators
        return new HighPrecision(left.num + right.num);
    }

    public static HighPrecision operator -(HighPrecision left, HighPrecision right)
    {
        // when we know the denominators are the same, can just subtract the numerators
        return new HighPrecision(left.num - right.num);
    }

    public override string ToString()
    {
        StringBuilder sb = new StringBuilder();

        // pull out the integer part
        BigInteger remain;
        BigInteger quotient = BigInteger.DivRem(num, HighPrecision.denom, out remain);
        int integerDigits = quotient.IsZero ? 1 : (int)BigInteger.Log10(quotient) + 1;
        sb.Append(quotient.ToString());

        int i = 0;
        BigInteger smallDenom = HighPrecision.denom/10;
        BigInteger tempRemain;

        // pull out all of the 0s after the decimal point
        while (i++ < HighPrecision.Precision && (quotient = BigInteger.DivRem(remain, smallDenom, out tempRemain)).IsZero)
        {
            smallDenom /= 10;
            remain = tempRemain;
            sb.Append('0');
        }

        // append the rest of the remainder
        sb.Append(remain.ToString());

        // truncate and insert the decimal point
        return sb.ToString().Remove(integerDigits + HighPrecision.Precision).Insert(integerDigits, ".");
    }
}

By far the most complicated part is the ToString() method, the rest is basic fraction operations with a shortcut or two thrown in.  Once we have the high precision number class, it is only a matter of plugging it in to the code for arctan:

public static HighPrecision GetPi(int digits)
{
    HighPrecision.Precision = digits;
    HighPrecision first = 4 * Atan(5);
    HighPrecision second = Atan(239);
    return 4 * (first - second);
}

public static HighPrecision Atan(int denominator)
{
    HighPrecision result = new HighPrecision(1, denominator);
    int xSquared = denominator * denominator;

    int divisor = 1;
    HighPrecision term = result;

    while (!term.IsZero)
    {
        divisor += 2;
        term /= xSquared;
        result -= term / divisor;

        divisor += 2;
        term /= xSquared;
        result += term / divisor;
    }

    return result;
}

Using this code, my laptop (2.5 GHz Intel Core i5) can churn out \(\pi\) to 1 million digits after the decimal point in a bit under 53 minutes. Producing the string representation takes another minute and forty seconds.  Using a fully homegrown solution (no BigInteger) with no attempt at optimization resulted in a runtime of about 8 hours, which is still in the realm of “reasonable” if you ask me.

The exercise I leave to the reader is to parallelize the above code for Atan.  After parallelizing with TPL, runtime on my work desktop (4 core, 2.8 GHz Intel Core i7) came down to about 17 minutes for 1 million digits.

Having done the calculation myself, I feel I have now earned the right to spend a lot of money on invest in this poster.  Or perhaps I will print my own.

Finally, I want to mention that this is still a very basic/nonoptimized approach compared to what leading researchers are doing.  Search around or check out this link for more gory details on how the pros do it.