Calculating PI in .NET

December 30, 2011

5 comments

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:

arctan(1)*4 as PI
  1. private static BigInteger CalculatePiSlow(int iterations, int digits) {
  2.     var mag = BigInteger.Pow(10, digits);
  3.     var sum1 = Enumerable.Range(0, iterations / 2).Select(i => i * 4 + 1).Aggregate(BigInteger.Zero, (sum, i) => sum += mag / i);
  4.     var sum2 = Enumerable.Range(0, iterations / 2).Select(i => i * 4 + 3).Aggregate(BigInteger.Zero, (sum, i) => sum += mag / i);
  5.     return 4 * (sum1sum2);
  6. }

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:

image

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:

image

No significant change. That same digit is still wrong. How about 10000 iterations?

image

This is better. The fourth digit is correct, but the fifth is not. Let’s get crazy and go for a million iterations:

image

Better still, but we’re not even at Math.PI. Ok, let’s go super-crazy and try 100 million iterations:

image

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:

image

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:

Arc Tangent for BigInteger
  1. static BigInteger ArcTan(int oneOverX, int iterations, int digits) {
  2.     var mag = BigInteger.Pow(10, digits);
  3.     var sum1 = Enumerable.Range(0, iterations / 2).Select(i => i * 4 + 1).Aggregate(BigInteger.Zero, (sum, i) => sum += mag / (i * BigInteger.Pow(oneOverX, i)));
  4.     var sum2 = Enumerable.Range(0, iterations / 2).Select(i => i * 4 + 3).Aggregate(BigInteger.Zero, (sum, i) => sum += mag / (i * BigInteger.Pow(oneOverX, i)));
  5.     return sum1sum2;
  6. }

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:

Faster PI calculation
  1. private static BigInteger CalculatePiQuick(int n, int digits) {
  2.     return 4 * (4 * ArcTan(5, n, digits) – ArcTan(239, n, digits));
  3. }

Running this with 20 iterations and 20 digits yields:

image

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):

image

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…

Add comment
facebook linkedin twitter email

Leave a Reply

Your email address will not be published.

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=""> <s> <strike> <strong>

*

5 comments

  1. Boris KozorovitzkyFebruary 3, 2012 ב 22:21

    Hi Pavel,

    First I would like to say that I enjoy reading your blog. It has many interesting and educating posts.

    I especially like this post because I also like playing with numbers.
    As a C++ and C# expert I would love to hear your comment on my (quite old) blog post here:
    http://alookonthecode.blogspot.com/2011/07/jitcompiler-optimization-of-fibonacci.html

    Where I try to calculate Fibonacci numbers.

    Thanks,
    Boris.

    Reply
  2. pavelyFebruary 4, 2012 ב 18:23

    Hi Boris,

    I’m glad you enjoy reading – I certainly enjoy writing!
    I’ve posted a comment on your blog post…

    Reply
  3. EvilGoatJune 6, 2013 ב 15:26

    Well, I think I see a pattern develop here as iterations need to be larger than the number of digits we look for. However I would like to know if there is a rough approximation how much larger or in fact how much it is in comparison to the digits (like for example digits*7/6 + 15) or something! 🙂

    GREAT POST BY THE WAY! 😉

    Reply
  4. Ross PresserJuly 5, 2016 ב 22:02

    This code can be improved by doing both Pow calculations in one lambda and by using PLINQ:

    {
    var mag = BigInteger.Pow(10, digits);
    var sum1 = Enumerable.Range(0, iterations / 2)
    .AsParallel()
    .Select(i => i * 4 + 1)
    .Select(i => mag / (i * BigInteger.Pow(oneOverX, i))
    - mag / ((i + 2) * BigInteger.Pow(oneOverX, i + 2))
    )
    .Aggregate(BigInteger.Zero,
    (sum, i) => sum += i)
    ;

    return sum1;
    }

    This got 4000 digits in 30.4 sec, on my Core i7 3.1Ghz.

    Reply
    1. Ross PresserJuly 5, 2016 ב 22:03

      excuse me, that was 10k digits, not 4k digits.

      Reply