Calculating PI, with computer or pencil

Computers have made the job of calculating ever more useless digits much easier - the current record is approximately 6.4 billion (you can download a very fast Windows based program for calculating millions of digits). Update: as of December 28, 2013, the record is 12.1 trillion digits of pi.

In the course of completing my high precision number class for Fractal eXtreme I had to create high precision versions of functions like exp, log, sin, cos, atan, sqrt, etc. I wanted to test these, and calculating PI to a few thousand digits seemed like a good way, so I scanned around for an appropriate formula for PI.

A little trigonometric thought will show that, since the tangent of a forty five degree angle is 1 (since sin and cos are equal at that point) then atan(1) is equal to forty five degrees, or PI/4 radians.

No problem, a naive implementor might think. I'll just calculate atan(1), perhaps using the handy Taylor series. The Taylor series for atan(x) is:

atan(x) = x - x^3/3 + x^5/5 - x^7/7 + x^9/9...*

therefore, simplifying for the case x equals 1:

PI/4 = atan(1) = 1 - 1/3 + 1/5 - 1/7 + 1/9...

Trouble is, the Taylor series for atan converges very slowly for values of x near one, and it converges extremely slowly when x is equal to one. My guess would be that if you wanted PI accurate to ten digits using this method, you would probably have to sum up at least a few billion terms. I would guess that this sequence is not terribly stable, so you would have to do all your math to quite high precision, just to produce a very small number of digits.

Not a good method.

A much better method is available however, courtesy of John Machin. It turns out that

PI/4 = atan(1) = 4 * atan(1/5) - atan(1/239)**

This is much more promising, because the Taylor series for atan(x) actually converges at a reasonable rate when x is equal to 1/5 or 1/239. Since each successive term is reduced by a factor of x squared, that means that with each term of atan(1/5) we gain about 1.4 digits. atan(1/239) converges much faster, with each term gaining about 4.75 digits. These swiftly diminishing terms have two advantages. Number one is that we need to do a relatively small number of terms to calculate PI to thousands of digits. Number two is that the calculations are quite stable. If a particular term is, say, 1e-50, then we know that all subsequent terms, and the sum of all subsequent terms, will be less than 1e-50. Therefore we know that at that point we have calculated atan to fifty digits of accuracy.

If we have a C++ high precision number class called InfPrec then a generic method for calculating atan(x) might look like this:

InfPrec atan(const InfPrec& x)
{
    InfPrec Result = x;
    InfPrec XSquared = x * x;
    
    InfPrec Term = x;
    int Divisor = 1;
    while (Term != 0)
    {
        Divisor += 2;
        Term *= XSquared;
        Result -= Term / Divisor;

        Divisor += 2;
        Term *= XSquared;
        Result += Term / Divisor;
    }
    return Result;
}

Unseen, but equally important, we need a way to initialize x with the values 1/5 and 1/239. This is probably best done by initializing the high precision numbers to one, and then dividing by an integer. We also probably want a way to print the results, preferably in decimal. This requires having a way to multiply a high precision number by an integer, so that we can repetitively multiply by ten and strip off the integer portion for printing.

So, calculating atan, and by extension PI, isn't too difficult at all.

However this method is both harder and slower than it needs to be. The hardest and slowest part of our algorithm is muliplying two high precision numbers together. This is mildly messy to write, especially in a high level language, and (unless you employ very sophisticated algorithms***) it is an O(n^2) algorithm. O(n^2) is just a fancy way of saying that doubling the number of digits will quadruple your calculation time.

We are using a generic atan routine. This is a handy thing to have, but overkill for this situation. We know that x will always be either 1/5 or 1/239. Therefore, instead of multiplying by a high precision number, we could simply divide by an integer - in this case, dividing by either 5^2 or 239^2. Normally I would not advocate replacing multiplication with division - it goes against everything we are taught in optimization school. But when the multiplication is going to involve millions of operations, and the division is going to involve just a few thousand, it can be a huge win.

Here's a special version of atan that handles the case that we need much more efficiently because XSquared is now an integer:

InfPrec ataninvint(int x)
{
    InfPrec Result = InfPrec(1) / x;
    int XSquared = x * x;

    InfPrec Term = Result;
    int Divisor = 1;
    while (Term != 0)
    {
        Divisor += 2;
        Term /= XSquared;
        Result -= Term / Divisor;
    
        Divisor += 2;
        Term /= XSquared;
        Result += Term / Divisor;
    }
    return Result;
}

With careful coding of these routines, a desktop computer can now calculate PI accurate to millions of digits! The code to calculate it is simply:

InfPrec PI = (ataninvint(5) * 4 - ataninvint(239)) * 4;

Pretty darned simple. On a 233Mhz Pentium II I was able to calculate one million digits of PI in just under six hours. More sophisticated algorithms can do this much faster but Machin's algorithm has a nice mixture of simplicity and efficiency.

What I think is really cool about this is how easily this algorithm works for manual calculation of PI. With a few sheets of paper, some patience, and a vague memory of long division, you can calculate PI to quite a high degree of accuracy. I estimate that manually calculating PI to thirty six digits would probably take a couple of days. Less than three hundred years ago****, that would have been enough for a world record, beating the efforts of a number of people who devoted much more time for worse results. They weren't using such efficient methods, because they hadn't yet been discovered.

When calculating PI by hand it makes sense to pre-multiply the two atan calls. This allows us to read off digits more easily, and makes the convergence more obvious. So, our formula is:

PI = 4 * (4 * atan(1/5) - atan(1/239))

or:

PI = 16 * atan(1/5) - 4 * atan(1/239)

Plugging this into our atan formula we get:

16 * atan(1/5) = 16/5 - 16/3*5^3 + 16/5*5^5 - 16/7*5^7 + 16/9*5^9...

and

4 * atan(1/239) = 4/239 - 4/3*239^3 + 16/5*239^5 - 16/7*239^7 + 16/9*239^9...

Calculating 16/5^n is particularly easy for people accustomed to dealing with powers of two. Since 1/5 equals 2/10, (1/5)^n can be rewritten as 2^n/10^n. In other words, a power of two, shifted right n digits. Multiplying by the additional sixteen makes our final calculation equivalent to 2^(n+4)/10^n.

Calculating the terms of atan(1/239) is a bit harder, because dividing by 239^2 is a bit messy. But, you need comparatively much fewer terms of atan(1/239) than of atan(1/5).

In the following tables, underlined digits represent digits which repeat. Noticing these repeating digits simplify the calculations tremendously because they frequently allow you to get an arbitrary number of digits of accuracy for a fixed amount of work.

16*atan(1/5)

4*atan(1/239)

Since addition is easier than subtraction, we're going to try to do as little subtraction as possible. Rather than subtracting the terms of atan(1/239), we simply negate them all so that we can then add them. Our answer is then the sum of all of the terms.

Half of those terms are still negative, so we won't be able to avoid all subtractions, but if we sum up all of the positive terms together, and all of the negative terms together, then we just have to subtract the sum of the negatives from the sum of the positives.

+3.200000000000000
+0.001024000000000
+0.000000910222222
+0.000000001008246
+0.000000000001233
+0.000000097666367
------------------
+3.201025008898068

-0.042666666666666
-0.000029257142857
-0.000000029789090
-0.000000000034952
-0.000000000000044
-0.016736401673640
-0.000000000001025
------------------
-0.059432355308274

PI = 3.201025008898068 - 0.059432355308274 = 3.141592653589794

The last digit should actually be a three rather than a four, but we have still managed to calculate PI to 14 decimal places!

You can download (27K) a sample Win32 program, with full, portable source code and generate a few thousand digits of PI yourself.

Thanks for reading. I hope you found this interesting. While you're here, take a look at Fractal eXtreme, which makes use of a much more heavily optimized and extended version of this InfPrec class. This allows fast and easy zooming in on many different types of fractals to virtually infinite magnification.

You can look here for more information on historical methods for calculating PI.

Drop me a line if you have any comments.

.Bruce Dawson.

Last modified: October, 2014

* Black magic #1. You'll just have to take it on faith that this series does actually sum up to atan(x). I do. Details and derivations can be found in searches for Taylor series or in calculus textbooks.

** Black magic #2. This formula (atan(1) equals 4 * atan(1/5) - atan(1/239)) is called the Machin Formula. You can read more about it, including how to easily derive it, on this page.

*** Although it is intuitively obvious that multiplication of high precision numbers must take n squared operations, it turns out that there are algorithms that improve the performance to take a number of operations proportional to n. There is considerable complication and overhead, but on numbers with millions of digits this can make multiplication run much much faster than the obvious method. Luckily we don't need this technique. For more information, read Knuth's excellent Art of Computer Programming.

**** Van Ceulen calculated PI to 35 digits c. 1600, apparently using Archimedes' method of inscribed and circumsribed polygons. In 1699 Sharp used PI/6 equals atan(30) to calculate 71 digits. Two years later Machin discovered and used the method described here and calculated 100 digits.

In the 1800's William Shanks used Machin's method to calculate 707 digits of PI, but only 527 of them were correct.

See this page for discussions of hand-calculations of pi.

Updated Dec 1, 2014, to remove an unused variable from ataninvint and to upload a new .zip file with the code updated to compile with modern C++ (iostream instead of iostream.h). I compiled it with VC++ 2013, with this command line:

cl SIMPLEPI.CPP /I"c:\Program Files (x86)\Microsoft Visual Studio 12.0\VC\include" /EHsc  -W4