In my previous post, I used a couple of algorithms to calculate Pi. With my implementation of the algorithms, to get to a precision of just 6 decimal places took well over a million iterations and over 3 seconds.

Clearly I was doing something wrong.

A few Google searches later I found this implementation by Julian M Bucknall. In this implementation, Julian creates a BigNumber class to handle the fact that, before .NET 4, there was no way to handle numbers to many digits of precision. His class uses fixed point arithmetic for storing data and performing calculations – effectively it works in base 4,294,967,296 (which is 2^{32}).

As .NET does now contain a BigInteger class (in .NET 4’s System.Numerics), I started searching for C# implementations that used this class. I couldn’t find anything in C#, but I did find a Java implementation that used a BigInteger class.

So, I decided to port it to C#, amending the BigInteger usages to match the .NET implementation syntax.

The formula used is the one that John Machin devised in 1706:

So, here is the .NET port of the Java implementation – it includes a method called InverseTan to be used in the above formula, which I won’t even try to understand at the moment. I’m just doing the port!

public class PiJavaPort { public static BigInteger InverseTan(int denominator, int numberOfDigitsRequired) { int demonimatorSquared = denominator*denominator; int degreeNeeded = GetDegreeOfPrecisionNeeded(demonimatorSquared, numberOfDigitsRequired); BigInteger tenToNumberPowerOfDigitsRequired = GetTenToPowerOfNumberOfDigitsRequired(numberOfDigitsRequired); int c = 2*degreeNeeded + 1; BigInteger s = BigInteger.Divide(tenToNumberPowerOfDigitsRequired, new BigInteger(c)); // s = (10^N)/c for (int i = 0; i < degreeNeeded; i++) { c = c - 2; var temp1 = BigInteger.Divide(tenToNumberPowerOfDigitsRequired, new BigInteger(c)); var temp2 = BigInteger.Divide(s, new BigInteger(demonimatorSquared)); s = BigInteger.Subtract(temp1, temp2); } Console.WriteLine("Number of iterations=" + degreeNeeded); // return s/denominator, which is integer part of 10^numberOfDigitsRequired times arctan(1/k) return BigInteger.Divide(s, new BigInteger(denominator)); } private static int GetDegreeOfPrecisionNeeded(int demonimatorSquared, int numberOfDigitsRequired) { //the degree of the Taylor polynomial needed to achieve numberOfDigitsRequired //digit accuracy of arctan(1/denominator). int degreeNeeded = 0; while ((Math.Log(2*degreeNeeded + 3) + (degreeNeeded + 1)*Math.Log10(demonimatorSquared)) <= numberOfDigitsRequired*Math.Log(10)) { degreeNeeded++; } return degreeNeeded; } private static BigInteger GetTenToPowerOfNumberOfDigitsRequired(int numberOfDigitsRequired) { var tenToNumberOfDigitsRequired = new BigInteger(1); // The following loop computes 10^numberOfDigitsRequired for (var i = 0; i < numberOfDigitsRequired; i++) { tenToNumberOfDigitsRequired = BigInteger.Multiply(tenToNumberOfDigitsRequired, new BigInteger(10)); } return tenToNumberOfDigitsRequired; } }

Finally, we need a method to put it all together – we need to implement the Machin formula:

public static string Calculate(int numberOfDigitsRequired) { numberOfDigitsRequired += 8; // To be safe, compute 8 extra digits, to be dropped at end. The 8 is arbitrary var a = BigInteger.Multiply(InverseTan(5, numberOfDigitsRequired), new BigInteger(16)); //16 x arctan(1/5) var b = BigInteger.Multiply(InverseTan(239, numberOfDigitsRequired), new BigInteger(4)); //4 x arctan(1/239) BigInteger pi = BigInteger.Subtract(a, b); var piAsString = BigInteger.Divide(pi, new BigInteger(100000000)).ToString(); var piFormatted = piAsString[0]+"."+piAsString.Substring(1,numberOfDigitsRequired-8); return piFormatted; }

Now to test the two implementations for performance – my Java port above vs. Julian M Bucknall implementation using his own BigNumber class.

First, my Java Port using .NET’s numerics:

Precision = 1000 digits Number of iterations needed = 3659 Elapsed Milliseconds = 96

Not bad! Much better than 1 million plus iterations for 6 decimal places that my pathetic attempt took!

Now to Julian M Bucknall’s implementation (using his own BigNumber class):

Precision = 1000 digits Number of iterations needed = 1603 Elapsed Milliseconds = 77

We have a winner.

Julian’s BigNumber class is optimised for this calculation – it can only divide a BigNumber by an integer. You cannot divide a BigNumber by another BigNumber. He states that if he had to implement such functionality he “would have discretely dropped the entire project”!

BigInteger does not have this limitation (and this functionality is used in the port).