I always loved mathematics. Although I’m certainly not a mathematician by profession, I’m always intrigued and inspired by math’s pureness and cleverness. One of the simplest and fascinating aspects of math is the number PI. Described simply as the ratio of a circle’s circumference to its diameter, it’s a constant with infinite digits after the decimal point and most importantly, non repeating (at least as far as I know).
There are many ways to calculate PI, as evident within the PI Wikipedia link. I wanted to see how I can get a large number of digits of PI with a C# program.
As every high school kid knows, PI in radians is equal to 180 degrees, or half a circle. And since the tangent of 45 (PI/4) degrees is 1, PI can be calculated like so:
PI/4 = arctan(1) or PI = 4 * arctan(1)
This looks simple enough. Just need to calculate the arc (inverse) tangent of 1. Of course, a method like Math.ATan is out of the question. It (and its kind) work with doubles, which gives a fixed 15 digit precision, while we want arbitrary precision.
How should we calculate the arc tangent of 1? One option is to use the Taylor series expansion (valid for |x|<= 1):
arctan(x) = x – x^3/3 + x^5/5 – x^7/7 + x^9/9 –…
The ^ sign indicates “to the power of”. If x is 1 (as is in our case), the series is simplified:
arctan(1) = 1 – 1/3 + 1/5 – 1/7 + 1/9 –…
Nothing to it, right?
But how do we calculate that? Again, using just doubles will get us nowhere better than Math.PI. We need an arbitrary sized floating point number. This does not exist in .NET (at least not yet). So, we’ll take the next best thing: the BigInteger type introduced in .NET 4. This is an arbitrarily large integer; all we have to do is use the above formula, but multiply everything by some large BigInteger, such as 10 to the 100th power, and interpret the final result as being larger by that factor. Let’s write a PI calculating method based on this series:
If you’re not familiar with the LINQ Aggregate operator, here’s a brief explanation: it’s the generic way operators like Sum and Average work: we use an “accumulator”, starting at BigInteger.Zero, and decide with the second argument (a delegate) what to do with the incoming number. The Sum operator does not work here because there is no overload that takes a BigInteger (although it would be simple enough to add such an extension method. I’ll leave it as an exercise for the interested reader).
iterations indicates the number of terms we wish to take into account and digits indicates the maximum digits we wish to calculate. Calling this method with 1000 iterations and just 20 digits yields the following:
I’ve printed Math.PI as a basic comparison. The result is disappointing – the fourth digit is wrong (“0” instead of “1”). Let’s bump up the iterations to 2000:
No significant change. That same digit is still wrong. How about 10000 iterations?
This is better. The fourth digit is correct, but the fifth is not. Let’s get crazy and go for a million iterations:
Better still, but we’re not even at Math.PI. Ok, let’s go super-crazy and try 100 million iterations:
This actually took a fair amount of time, but the result is still disappointing. Clearly, something must change.
The problem with this series is that it converges slowly for large x (close to 1). In fact, for very small x values, it converges pretty fast. We need a new formula.
In 1706, John Machin found this formula to calculate PI:
At first, this may seem worse: now we have to calculate two arctangents! But, as the numbers are much smaller than one, we should get a better rate of convergence.
First, we’ll have to write a general arc tangent function, as none exists in BigInteger. This is a simple variant of the above calculation that needs to raise to the appropriate power – a simple enough job as BigInteger natively supports it:
The last part uses the BigInteger.Pow method to raise to the appropriate power. Now we can write the new PI calculation method like so:
Running this with 20 iterations and 20 digits yields:
The result already surpasses Math.PI, and is calculated in a fraction of a second. We can try 1000 for the iterations and digits and get the first 1000 digits (or so) of PI pretty quickly (less than half a second without any optimizations):
The last 6 digits are in fact wrong (we need a few more iterations to get them right).
Now here’s a challenge for you: memorize as many digits as you can…