Kash Farooq's software development blog

.NET Developer

Posts Tagged ‘Pi’

Calculating Pi in C# part 3 – using the BigRational class

Posted by Kash Farooq on August 1, 2011

In an earlier post I ported a Java implementation of a Pi calculator – my port used the BigInteger class that lives in the System.Numerics assembly of .NET 4. I also performance tested my port against an implementation by Julian M Bucknall that used his own BigNumber class. Julian’s implementation won that performance test, coming in at an impressive 70-80 milliseconds to calculate Pi to 1000 decimal places.

Now it is time for me to try an implementation alone!

I’m going to use the following John Machin formula:

John Machin's Pi formulaArctan can be calculated using the following Taylor series:

ArcTan calculated using the Taylor SeriesRather than creating arbitrarily big numbers with BigInteger, I need to create numbers with arbitrarily large precision.   I need a BigRational. And one exists, but it did not make it into the .NET 4 framework. You can download BigRational form CodePlex.

So, here is my implementation of arctan:

public class ArcTanWithBigRational {
    public ArcTanWithBigRational() {
        Iterations = 0;
    }

    public int Iterations;

    public BigRational Calculate(BigRational x, int maxNumberOfIterations, int precision) {
        bool doASubtract = true;
        BigRational runningTotal = x;
        int count = 0;
        var divisor = new BigInteger(3);
        while (count < maxNumberOfIterations) {
            BigRational current = BigRational.Pow(x, divisor);
            current = current/divisor;
            if (doASubtract) {
                runningTotal = runningTotal - current;
            }
            else {
                runningTotal = runningTotal + current;
            }
            doASubtract = !doASubtract;
            count++;
            divisor = divisor + 2;
            if (WeHaveEnoughPrecision(current, precision)) {
                Iterations = count;
                break;
            }
        }
        return runningTotal;
    }

    private static bool WeHaveEnoughPrecision(BigRational current, int precision) {
        return current.GetFractionPart().ToString().Length > precision+2; //extra 2 digits to ensure enough precision
    }
}

Now to perform the calculation:

public string Calculate(int numberOfDigits) {
    var arcTanA = new ArcTanWithBigRational();
    var arcTanB = new ArcTanWithBigRational();
    var a = 16 * arcTanA.Calculate(BigRational.Divide(1, 5), 1000, numberOfDigits);
    var b = 4 * arcTanB.Calculate(BigRational.Divide(1, 239), 1000, numberOfDigits);
    var pi = a - b;
    return BigRationalPiFormatter.Format(pi, numberOfDigits);
}

And the results:

Precision = 1000 digits
Number of iterations = 923
Elapsed Milliseconds = 1192

Not bad! No where near as fast as the other implementations in my previous post, but not too bad and with far fewer iterations.

I can try one final optimization to try and get my implementation running faster. I have a dual-core laptop so let’s try some Parallel Tasks:

public string Calculate(int numberOfDigits) {
    var a = new ArcTanWithBigRational();
    var b = new ArcTanWithBigRational();
    Task<BigRational> task1 = Task<BigRational>.Factory.StartNew(
                            () => a.Calculate(BigRational.Divide(1, 5),1000,numberOfDigits));
    Task<BigRational> task2 = Task<BigRational>.Factory.StartNew(
                            () => b.Calculate(BigRational.Divide(1, 239), 1000, numberOfDigits));

    var pi = 16 * task1.Result - 4 * task2.Result;

    return BigRationalPiFormatter.Format(pi, numberOfDigits);
}

The results this time:

Precision = 1000 digits
Number of iterations = 923
Elapsed Milliseconds = 1099

A little faster, but nothing to get excited about.

For completeness, here is my routine to format the BigRational Pi into a string:

public class BigRationalPiFormatter {
    public static string Format(BigRational pi, int numberOfDigits)
    {
        BigInteger numeratorShiftedToEnoughDigits =
                       (pi.Numerator * BigInteger.Pow(new BigInteger(10), numberOfDigits));
        var bigInteger = numeratorShiftedToEnoughDigits / pi.Denominator;
        string piToBeFormatted = bigInteger.ToString();
        var builder = new StringBuilder();
        builder.Append(piToBeFormatted[0]);
        builder.Append(".");
        builder.Append(piToBeFormatted.Substring(1, numberOfDigits - 1));
        return builder.ToString();
    }
}

Posted in .NET, Algorithms | Tagged: , , , | 2 Comments »

Calculating Pi in C# part 2 – using the .NET 4 BigInteger class

Posted by Kash Farooq on July 23, 2011

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 232).

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:

John Machin's Pi formulaSo, 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).

Posted in .NET, Algorithms | Tagged: , , , | Leave a Comment »

Calculating Pi in C# with series algorithms

Posted by Kash Farooq on July 17, 2011

People write applications to calculate Pi to thousands of decimal places – I decided to investigate how they do that.

Google took me to a post discussing using recursion to calculate Pi. Quite ironic considering the post is on Stack Overflow… Anyway, the post did provide me with plenty of links and algorithm names for me to refine my search and find a way to do it without recursion.

I soon found this: a collection of series algorithms to calculate Pi.

So, I picked an algorithm and started coding.

To get going quickly, I decided to just use the built in .NET value types (i.e. decimal, double). These clearly are not good enough if you want to calculate Pi to hundreds of decimal places, but they would do for now.

To check my algorithm implementation worked, I got Pi to a few thousand decimal places from the interweb :

public static class Pi {
    public const string PiAsString = "3.14159265358979323846264338327...etc";
    public static string Get(int numberOfDecimalPlaces) {
        return PiAsString.Substring(0, numberOfDecimalPlaces + 2);
    }
}

Leibniz-Gregory-Madhava Series

You can calculate Pi using the following:

In the code below, my base class property “PiCalculatedToRequiredPrecision” is, as the name suggests, there stop if I’ve reached my target number of digits. The “KeepGoing” property is used to bail out of the loop if I’ve done far too many iterations (which I keep track of in the Iteration property) and I have still not managed to get the required number of digits.

public class LeibnizGregoryMadhava:SeriesAlgorithmBase {
    public LeibnizGregoryMadhava(int numberOfDecimalPlaces) : base(numberOfDecimalPlaces) {}

    // Pi/4 = 1 - 1/3 + 1/5 - 1/7 + 1/9 - 1/11 .....

    public override decimal Calculate() {
        CurrentPi = 1;
        var multiplier = 1;
        var nextNumber = 3;

        decimal result = 1;
        while (KeepGoing) {
            multiplier = multiplier * -1;
            result = result + multiplier * (decimal)1 /(nextNumber);
            CurrentPi = result*4;
            if (PiCalculatedToRequiredPrecision)
            {
                break;
            }
            Iteration++;
            nextNumber += 2;
        }
        return CurrentPi;
    }
}

And now to run my algorithm a few times specifying a different precision each time:

LeibnizGregoryMadhava took 117 iterations and 6 milliseconds to calculate PI to 2 decimal places
LeibnizGregoryMadhava took 1686 iterations and 3 milliseconds to calculate PI to 3 decimal places
LeibnizGregoryMadhava took 10792 iterations and 22 milliseconds to calculate PI to 4 decimal places
LeibnizGregoryMadhava took 136119 iterations and 290 milliseconds to calculate PI to 5 decimal places
LeibnizGregoryMadhava took 1530010 iterations and 3358 milliseconds to calculate PI to 6 decimal places

Over 1.5 million iterations (and 3.4 seconds) to get to 6 decimal places!!

Time to try another algorithm.

Euler Series

public class Euler:SeriesAlgorithmBase {
    //pi^2 / 6 = 1 + 1/4 + 1/9 + ..... + 1/k^2

    public Euler(int numberOfDecimalPlaces) : base(numberOfDecimalPlaces) {}
    public override decimal Calculate() {
        decimal result = 0;
        Iteration = 1;
        while (KeepGoing)
        {
            long squared = Iteration * Iteration;
            result = result + (decimal)1/squared;
            CurrentPi = (decimal) Math.Sqrt((double) (result * 6));
            if (PiCalculatedToRequiredPrecision())
            {
                break;
            }
            Iteration++;
        }
        return CurrentPi;

    }
}

And the results:

Euler took 600 iterations and 74 milliseconds to calculate PI to 2 decimal places
Euler took 1611 iterations and 3 milliseconds to calculate PI to 3 decimal places
Euler took 10307 iterations and 22 milliseconds to calculate PI to 4 decimal places
Euler took 359863 iterations and 819 milliseconds to calculate PI to 5 decimal places
Euler took 1461054 iterations and 3319 milliseconds to calculate PI to 6 decimal places

Hmmm. Fewer iterations (just!) and 3.3 seconds to get to 6 decimal places.

I’m not going to break any world records with this code. Clearly I need to do more research.

Related posts

Posted in .NET, Algorithms | Tagged: , , , | 3 Comments »