Calculating Pi

One programming task that I had never taken on (up until this morning, anyway), was a pi calculator. Finding pi is not an overwhelmingly daunting task; one can get an extremely precise approximation from the System.Math.PI constant. However, that’s no fun. The real fun is in rolling your own pi calculator. As it turns out, this was no beginner project. It was math-heavy, theory-heavy and CPU-heavy. Oh, and it was a blast!

The trouble was that at first, I had no idea where to begin. In C#, there is no “arbitrary precision” decimal type; the decimal class can hold 28 or so significant figures (remember the calculator tutorial?). No self-respecting pi calculator would ever stop at 28 significant figures. 28,000 would be a bit better! The first challenge in coding the pi calculator would be coming up with a way to store an unlimited number of significant figures. I’d read about some pi calculator authors implementing “BigFraction” or “BigRational” classes based off of the fairly new BigInteger (arbitrarily-sized .NET integer) class, so I thought I’d take a crack at writing one of those. I also decided that before I wrote a single line of code, I needed to choose a pi calculation formula.

The requirements for a formula are simple: that it converges relatively quickly and requires relatively few calculations per iteration. Unfortunately, the use of a “BigRational” class adds another requirement: no fractional powers of numbers that have irrational fractional powers. That’s a real drag, because it means we can’t use the Chudnovsky Algorithm, which is the Holy Grail of fast pi calculation algorithms:

The trouble is with the 3k + 3/2 part; 640320 is not a perfect square (it is almost 800.2^2, but not quite).

The trouble is with the 3k + 3/2 part; 640320 is not a perfect square (it is almost 800.2^2, but not quite).

With the really good algorithm out of the way, I started trying to find another good algorithm. My first thought took me back to my trig identity days a few years ago (okay, the past three years of high school):

So then if you take the deriv... actually, no.

So then if you take the derivat… actually, no.

This is actually a sound, correct thought. In fact it has a name: the Gregory-Leibniz formula. This will eventually converge to pi. However, after following a link off the handy, dandy Wikipedia page I used to decide which algorithm I would use in my program, I learned that it takes five billion iterations of arctan(1) (which I will get to in a moment) to produce ten correct digits of pi. YIKES! It was a good thought.

Anyway, I eventually decided on the original Machin formula, from 1706:

Although it was created in 1706, Wikipedia says it is still one of the fastest-converging pi-calculation algorithms!

Despite its age, it converges quickly and is relatively easy to compute.

Awesome. We have a formula. Now we need to figure out how an arctangent can be computed. It is actually rather simple, especially with programmatic help:

See. Not that bad.

See. Not that bad.

Alright! We have the formula for computing pi and we have the formula for computing arctangent, so now we’re done. Not so fast. Remember when I mentioned the BigRational class? As it turns out, writing a BigRational class is the bulk of this project. I won’t bore you with too many tiny details of the way I implemented it; I will tell you that it was designed for simplicity over speed. It definitely needs some refining, but it works.

The basic concept behind the BigRational class is simple: it is a fraction. It has a numerator and a denominator. The numerator and denominator are arbitrarily-precise, so the entire ratio is arbitrarily-precise. The class has a bunch of methods to help the numerator and denominator do their fraction-y duty:

As it turns out, a lot of help.

As it turns out, a lot of help.

Most of the methods listed in the class diagram are pretty self-explanatory; I won’t go into detail here (but you can download the source code to the project at the bottom of the page). All the operators help the BigRationals play well together, and the * operator has an overload that allows “scalar” multiplication. That is, it allows multiplication by a BigInteger.

The two methods that I will delve into are the Reduce method and the ToDecimalString method. The Reduce method reduces a BigRational into its simplest form. It implements Euclid’s Algorithm, which is a means of reducing a fraction. It works in steps:

  1. The larger of the numerator and the denominator is set as the “big number.” The smaller is set as the “small number.”
  2. The program finds the remainder when the big number is divided by the small number.
  3. If the remainder is not zero, the small number becomes the big number and the remainder from step two becomes the small number.
  4. If the remainder is zero, both the numerator and denominator may be divided by the small number. At this point, the small number is the greatest common factor of the numerator and denominator.
  5. Steps 2, 3 and 4 are repeated until step four terminates with the small number equal to one. At this point, the fraction is in simplest form.

In code, the above steps look like this:

public void Reduce()
{
    // Set up big and small number
    BigInteger bigNumber, smallNumber = 0;
    if (Numerator > Denominator)
    {
        bigNumber = Numerator;
        smallNumber = Denominator;
    }
    else
    {
        bigNumber = Denominator;
        smallNumber = Numerator;
    }
 
    // Now divide a bunch of times
    while (smallNumber != 1)
    {
        // Find the remainder
        BigInteger rem = (bigNumber % smallNumber);
        if (rem != 0)
        {
            bigNumber = smallNumber; // Set up the numbers for the next iteration
            smallNumber = rem;
        }
        else
        {
            // We can divide both the numerator and denominator by the previous remainder
            Numerator /= smallNumber;
            Denominator /= smallNumber;
 
            // Re-assign bigNumber and smallNumber
            if (Numerator > Denominator)
            {
                bigNumber = Numerator;
                smallNumber = Denominator;
            }
            else
            {
                bigNumber = Denominator;
                smallNumber = Numerator;
            }
        }
    }
 
    // One last thing -- if the numerator is positive and the denominator is negative, the numerator needs to become negative and the denominator needs to become positive.
    // If they're both negative, they should both be positive
    if (Denominator < 0 && Numerator > 0 || Numerator < 0 && Denominator < 0)
    {
        Numerator = -Numerator;
        Denominator = -Denominator;
    }
}

The only change between the pseudo-code and actual implementation was the addition of the small if-statement at the bottom. I noticed that the denominator would sometimes be negative instead of the numerator being negative. That would probably throw off some calculations, so I made sure that only the numerator of the BigRational could ever be negative.

The ToDecimalString method converts a BigRational into its decimal equivalent — a handy helper when it comes to calculating pi (since it’s no fun unless it starts with 3.14159…). I considered two different implementation styles:

  1. Multiplying the entire fraction by 10^however many digits I wanted.
  2. Creating my own long division algorithm that would return a string of digits of a specified length.

Contrary to what seems sane, I chose the second method. I suppose it was partially because I had always wanted to write an actual long division algorithm, but by any means, it took some head-scratching at first. Long division isn’t something I have to do very much, so I started by doing out a few long division problems. As it turns out, there is a nice pattern. The non-decimal part of the number may be calculated in one fell swoop by using integer division. It may then be committed to a string. After that, a decimal point should be committed to the string. Then, a “new numerator” should be created. Its value is the remainder from the previous iteration multiplied by ten. The next decimal digit is equal to the “new numerator” minus the remainder when it is divided by the denominator of the fraction, all divided by the denominator of the fraction. This process can be repeated until the desired number of digits is reached.

Okay, it makes more sense in code. Wow.

public string ToDecimalString(BigInteger decimalDigits)
{
    StringBuilder rv = new StringBuilder();
    this.Reduce(); // Go as fast as possible
    BigInteger remainder = Numerator % Denominator;
 
    // Get the non-decimal part
    rv.Append(((Numerator - remainder) / Denominator).ToString() + ".");
    remainder = Numerator % Denominator;
 
    BigInteger newNum = remainder * 10;
 
    // Now get the decimal part
    for (BigInteger i = 0; i < decimalDigits; i++)
    {
        rv.Append(((newNum - (newNum % Denominator)) / Denominator).ToString()); // This literally just does long division
        newNum = (newNum % Denominator) * 10;
    }
 
    return rv.ToString();
}

The only real advantage this method presents is the lack of any string parsing. It is all basic use of a StringBuilder. Plus, it seems to work very well!

With the BigRational class done, the next step was the actual pi calculation. The first step was to implement an arctangent method, in accordance with the formula I mentioned earlier.

public static BigRational ArcTangent(BigRational input, BigInteger iterations)
{
    BigRational retVal = input;
    for (BigInteger i = 1; i < iterations; i++)
    {
        // arctan(x) = x - x^3/3 + x^5/5 ...
        // = summation(n->infinity) (-1)^(n) * x^(2n+1)/(2n+1)
        BigRational powRat = input.Pow((2 * i) + 1);
        retVal += new BigRational(powRat.Numerator * (BigInteger)Math.Pow(-1d, (double)(i)), ((2 * i) + 1) * powRat.Denominator);
        if (i % 100 == 0)
        {
            Console.WriteLine("ArcTangent {0}: {1}/{2} iterations complete.", input, i, iterations); // Status update.
        }
    }
 
    return retVal;
}

Wait. There’s a parameter called “iterations.” How do we figure out how many iterations are necessary? This page provided the answer. As it turns out…

I find it's a good idea to add twenty or so digits to the end to make sure everything is correct.

I find it’s a good idea to add twenty or so digits to the end to make sure everything is correct.

In other words…

public static BigInteger GetConversionIterations(BigInteger digits, BigRational q)
{
    return (BigInteger)((double)digits / (2 * Math.Log10((double)q.GetReciprocal())));
}

That was quick. That also finished laying the foundation for actually calculating pi. Phew. That was a lot of keystrokes.

public static BigRational GetPi(BigInteger numDigits)
{
    // pi = 16 arctan(1/5) − 4 arctan(1/239)
    BigRational oneFifth = new BigRational(1,5);
    BigRational oneTwoThirtyNine = new BigRational(1, 239);
    BigRational arcTanOneFifth = PiCalc.ArcTangent(oneFifth, PiCalc.GetConversionIterations(numDigits + 1, oneFifth)); // Start computing
    BigRational arcTanOneTwoThirtyNine = PiCalc.ArcTangent(oneTwoThirtyNine, PiCalc.GetConversionIterations(numDigits + 1, oneTwoThirtyNine));
 
    return (arcTanOneFifth * 16) - (arcTanOneTwoThirtyNine * 4);
}

This code block implements Machin’s formula as discussed above. An improved version of this program could take advantage of multi-core computing by calculating one arctangent on one thread and the other arctangent on another. More formulas should also be tested to determine which is the fastest; this program really slows down as the iteration count gets higher because the fractions get really, really, really huge. Of course, this program was never going to break any records anyway because it was written in C#. JIT-compiled, super high-level languages aren’t exactly speed demons. Then again, my Phenom II x4 (overclocked to 3.7 gHz) is also beginning to show its age. I’d be curious to see how performance compares on one of the latest-generation i7s.

Although it is sluggish, it does seem to compute pi properly. Here are the first 10,001 digits:

Although it looks like a lot of pi, it's not a lot of pi.

Although it looks like a lot of pi, it’s not a lot of pi.

So that’s it. It’s done. It works. Now it just needs to run for a few days so I can brag to my friends! The funny thing about calculating pi is how although we believe pi is completely irrational and will never end, we’ve calculated enough digits that it might as well end. If my calculations are correct, it only takes about 80 digits of pi to calculate the volume of the observable universe to the nearest cubic meter. Thus, ten thousand would calculate it to the nearest 1/10^9920 cubic meter, whatever that would be. That’s some mind-blowing precision. Imagine what a million digits would do.

Download this Project (MIT License)

Leave a Reply

Your email address will not be published. Required fields are marked *


5 − = one

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>