Estimating pi

Calculating Pi

package pi;

import java.util.*;
import java.math.BigInteger;
import java.math.BigDecimal;

/**
* Calculates pi to an arbitrary number od decimal places
* @author BotRejectsInc.
*/
public class Main
{

/** Creates a new instance of Main */
public Main()
{
}

public static void main(String[] args)
{
    int terms = 1;
    Scanner sc = new Scanner(System.in);

    while(true)
    {
        System.out.println("Enter a positive integer { 1, 2, 3, ...} or 0 to exit:");
        terms = sc.nextInt();
        if(terms < 1)
        { break; }
        BigDecimal pi = CalculatePi(terms);
        System.out.println(pi);
        //System.out.println("As an approximating double: " + pi.doubleValue());
    }

    sc.close();
}

private static BigDecimal CalculatePi(int terms)
{
    //Need to use BidDecimals for this!
    int decimalPlaces = 1000;
    BigDecimal inverse_pi, Pi, numerator, divisor;
    BigDecimal[] sigma = new BigDecimal[terms];
    BigDecimal[] term = new BigDecimal[terms];

    //This coeeficient multiplies the final sum of terms
    BigDecimal root_2 = Sqrt_2(200);
    BigDecimal coefficient = root_2.multiply(BigDecimal.valueOf(2)).divide(BigDecimal.valueOf(9801), 1000,
BigDecimal.ROUND_HALF_UP);

    //Calculate sum of terms
    for(int i = 0; i < terms; i++)
    {
        //Calculate next term
        try
        {
            numerator = new BigDecimal( LargeFactorial(4*i).multiply(BigInteger.valueOf(1103 + 26390*i)) );
            divisor = new BigDecimal( LargeFactorial(i).pow(4).multiply( BigInteger.valueOf(396).pow(4*i) ) );
            term[i] = ( numerator.divide(divisor, decimalPlaces, BigDecimal.ROUND_HALF_UP) );
        }
        catch( IllegalArgumentException e )
        {
            System.out.println(e.getMessage());
        }
        //Summate terms and store cumulative sum in the array sigma
        if(i == 0)
        {
            sigma[i] = term[i];
        }
        if(i > 0)
        {
            sigma[i] = term[i].add(sigma[i-1]);
        }
    }
    //Final algebra multiply sum by coefficient and invert
    inverse_pi = coefficient.multiply(sigma[terms-1]);
    BigDecimal unity = new BigDecimal("1");
    Pi = unity.divide(inverse_pi, decimalPlaces, BigDecimal.ROUND_HALF_UP);

    return Pi;
}

//Calculates factorials of large values using BigInteger
private static BigInteger LargeFactorial(int n)
throws IllegalArgumentException
{
    if(n == -1)
    {
        throw new IllegalArgumentException("Negative factorial not defined");
    }
    BigInteger factorial = new BigInteger("1");
    if(n == 0 || n == 1)
    {
        return factorial;
    }
    while(n > 1)
    {
        BigInteger N = BigInteger.valueOf(n);
        factorial = factorial.multiply(N);
        //System.out.println(factorial);
        n--;
    }
    return factorial;
}

private static BigDecimal Sqrt_2(int terms)
{
    BigDecimal root_2;
    BigDecimal[] buffer = new BigDecimal[terms+1];
    BigDecimal initial = new BigDecimal("1.4");
    buffer[0] = initial;
    for(int i = 1; i <= terms; i++)
    {
        buffer[i] = buffer[i-1].divide(BigDecimal.valueOf(2), terms*5,
BigDecimal.ROUND_HALF_UP).add(BigDecimal.valueOf(1).divide(buffer[i-1], terms*5,
BigDecimal.ROUND_HALF_UP));
    }
    root_2 = buffer[terms];
    return root_2;
}

}

Pi is an irrational number, meaning that it apparently has an infinite number of decimal places. My pocket
calculator gives pi as 3.141592654, a reasonable approximation! As someone who programs computers for a
hobby, I like to give myself little programming challenges and tasks to do. The aim of this program was to
calculate pi to at least 100 decimal places (d.p.) accuracy (essentially to an arbitrary accuracy, within the limits of
computation power on a desktop PC). (This kept me busy for a couple of afternoons and was very interesting
and educational!) The solution I obtained demonstrates the use of the Java BigInteger and BigDecimal classes,
the use of iteration and summation of terms in mathematical series. For this purpose I used the following
approximation formula for pi (the Ramanujan formula):

Note that this is a series approximation - we calculate a series of terms for k = 0, k = 1, k = 2, ..., etc. and
then add the terms together. The more decimal places of accuracy we require, the more terms we must
summate. This series is a good approximation for computation, since it converges very rapidly on the true
value of pi. We may wish to add tens or even hundreds of these terms together.

This looks simple enough, until you realise the problem caused by the factorials (! is factorial or 'shriek').
Recall that n! = n x n - 1, x, n-2, ..., x 2, so for example, 5! = 5 x 4 x 3 x 2 = 120 and 10! = 3 628 800 (0! is
set equal to 1 by definition). These numbers rapidly get very big! The highest my pocket calculator can
go is 69! = 1.711224524 E+98. For the 20th term, k = 20, we already need to compute 80!

It would be straightforward to multiply Java (or C#) integers (int) together to compute n! except that an int
is 32 bit and so can only hold a positive number as large as 2^(32)-1, that is 4294967295. A long is a 64
bit representation of whole numbers and so can go much higher, up to about 1.845 E+19, but is still not
large enough! Using these data types can result in overflow errors: if the number becomes to large for the
int or long container then only part of the number is stored and it may have a negative value. To avoid this
Java and C# come equipped with a BigInteger class which can store and manipulate whole numbers of
arbitrary length. Using this class, we can write a function to compute arbitrarily large factorials, such as:

//Calculates factorials of large values using BigInteger
private static BigInteger LargeFactorial(int n)
throws IllegalArgumentException
{
   if(n == -1)
   {
       throw new IllegalArgumentException("Negative factorial not defined");
   }
   BigInteger factorial = new BigInteger("1");
   if(n == 0 || n == 1)
   {
       return factorial;
   }
   while(n > 1)
   {
       BigInteger N = BigInteger.valueOf(n);
       factorial = factorial.multiply(N);
       //System.out.println(factorial);
       n--;
   }
   return factorial;
}

Generally speaking, the factorial of a negative number is not defined (though it can be in advanced
mathematics) and so we need to ensure that only zero or a positive integer is passed to the function.
There are a number of ways to do this, I chose to have it throw an IllegalArgumentException if it is passed
a negative argument.

Testing this function it certainly agrees with my calculator up to 69! and so presumably is accurate for
larger factorials. Here are a few test results:

The factorial of 10 is:   3628800

The factorial of 12 is:   479001600

The factorial of 13 is:   6227020800

The factorial of 14 is:   87178291200

The factorial of 15 is:   1307674368000

The factorial of 16 is:   20922789888000

The factorial of 17 is:   355687428096000

The factorial of 18 is:   6402373705728000

The factorial of 20 is:   2432902008176640000

The factorial of 30 is:   265252859812191058636308480000000

The factorial of 40 is:   815915283247897734345611269596115894272000000000

The factorial of 50 is:   30414093201713378043612608166064768844377641568960512000000000000

The factorial of 69 is:
171122452428141311372468338881272839092270544893520369393648040923257279754140647424000000000000000

The factorial of 100 is:
93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

The factorial of 200 is:
78865786736479050355236321393218506229513597768717326329474253324435944996340334292030428401198462390417721213891963883025764279024263710506192662495282993111346285727076331723739698894392244562145166424025403329186413122742829485327752424240757390324032125740557956866022603190417032406235170085879617892222278962370389737472000000
0000000000000000000000000000000000000000000

So we can calculate the factorials we need. What about the square-root of 2 which appears in the
coefficient that multiplies the sum of terms? This is an irrational number, like pi, with an apparently infinite
number of d.p. Using standard elementary data types, such as doubles or floats, again is insufficient - the
values overflow. Using doubles introduces an error into the 15th d.p. so we can not use something like:

double a = (2*Math.sqrt(2)) / 9801;

Java has another class for fractions, the BigDecimal class, but this does not support a square-root
method! The solution I used was to approximate root 2 using the following iterative formula:

This is an iterative formula that needs an initial 'guess' to get going. The value of the initial guess is not
critical, the iteration converges to root 2, but using a closer value, like 1, will mean that it converges more
quickly. I used 1.4 as an initial value. This formula does indeed converge very rapidly. Again, it is an
approximation, so we need to make sure that it is at least as accurate as we need to estimate pi to our
required degree of accuracy. I wrote a function which takes the number of terms or iterations as an
argument:

private static BigDecimal Sqrt_2(int terms)
{
   BigDecimal root_2;
   BigDecimal[] buffer = new BigDecimal[terms+1];
   BigDecimal initial = new BigDecimal("1.4");
   buffer[0] = initial;
   for(int i = 1; i <= terms; i++)
   {
       buffer[i] = buffer[i-1].divide(BigDecimal.valueOf(2), terms*5,
BigDecimal.ROUND_HALF_UP).add(BigDecimal.valueOf(1).divide(buffer[i-1], terms*5,
BigDecimal.ROUND_HALF_UP));
   }
   root_2 = buffer[terms];
   return root_2;
}

The result of each iteration is stored as a BigDecimal in the array buffer[], and these are then summated.
This was done as the BigDecimal, like the BigInteger type, is immutable (though references can be
reassigned as in the LargeFactorial() function above which did not use an array to store intermediate
values). It also means we can inspect any of the terms in future should we wish to do so. Be careful when
carrying out division with the BigDecimal class to specify a limit to the number of decimal places of
accuracy, as the division function can not cope with an endless string of digits (again ensure that the limit
used is at least sufficient for our requirements, by setting it as required in this case). I am also following
standard procedure of rounding halfs up.

We can test the results against the official value of root 2:

1.41421356237309504880168872420969807856967187537694807317667973799 ...

I obtained:

1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415728

with 100 iterations. This seems to be working fine, so lets test the final program, which is included below
as a single listing (you can run this as a single class and file console app in an IDE like NetBeans, hence
the use of static functions).

The official value of pi to 100 d.p. is:

3.1415926535 8979323846 2643383279 5028841971 6939937510 5820974944 5923078164
0628620899 8628034825 3421170679 8214808651 3282306647 0938446095 5058223172 5359408128
4811174502 8410270193 8521105559 6446229489 5493038196 4428810975 6659334461 2847564823
3786783165 2712019091 4564856692 3460348610 4543266482 1339360726 0249141273 7245870066
0631558817 4881520920 9628292540 9171536436 7892590360 0113305305 4882046652 1384146951
9415116094 3305727036 5759591953 0921861173 8193261179 3105118548 0744623799 6274956735
1885752724 8912279381 8301194912 9833673362 4406566430 8602139494 6395224737 1907021798
6094370277 0539217176 2931767523 8467481846 7669405132 0005681271 4526356082 7785771342
7577896091 7363717872 1468440901 2249534301 4654958537 1050792279 6892589235 4201995611
2129021960 8640344181 5981362977 4771309960 5187072113 4999999

With just 100 terms of the Ramanujan formula I obtained:

100 terms:

3.1415926535 8979323846 2643383279 5028841971 6939937510 5820974944 5923078164
0628620899 8628034825 3421170679 8214808651 3282306647 0938446095 5058223172 5359408128
4811174502 8410270193 8521105559 6446229489 5493038196 4428810975 6659334461 2847564823
3786783165 2712019091 4564856692 3460348610 4543266482 1339360726 0249141273 7245870066
0631558817 4881520920 9628292540 9171536436 7892590360 0113305305 4882046652 1384146951
9415116094 3305727036 5759591953 0921861173 8193261179 3105118548 0744623799 6274956735
1885752724 8912279381 8301194912 9833673362 4406566430 8602139494 6395224737 1907021798
6094370277 0539217176 2931767523 8467481846 7669405132 0005681271 4526356082 7785771342
7577896091 7363717872 1468440901 2249534301 4654958537 1050792279 6892589235 4201995611
2129021960 8640344181 5981362977 4771309960 5187072113 4999999

(8372978049 9510597317 3281609631 9125322977 4937741372 6512882026 8619471299
1299221042 1681419479 8705721777 9038338672 9847360285 5172958903 9759373423 2587356335
2043217855 3912536912 7618906019 4054262794 9728492904 0581146064 0009566664 5502282746
007 )

Accurate to 100 d.p. plus a few, and the calculation took only a fraction of a second!

With fewer terms less accuracy is obtained:

20 terms:

3.1415926535 8979323846 2643383279 5028841971 6939937510 5820974944 5923078164
0628620899 8628034825 3421170679 8214808651 3282306647 0938446095 5058223172 5359408128
481117450 (rest inaccurate)

30 terms:

3.1415926535 8979323846 2643383279 5028841971 6939937510 5820974944 5923078164
0628620899 8628034825 3421170679 8214808651 3282306647 0938446095 5058223172 5359408128
4811174502 8410270193 8521105559 6446229489 5493038196 4428810975 6659334461 2847564823
37867831 (rest inaccurate).

Developing the Program Further

The program could be modified by enabling the user to obtain pi to a requested accuracy: the user could
input the number of decimal places required. The solution here would involves looping indefinitely,
summating the terms in the series and checking the answer each time (with n terms) with the previous
answer (with n - 1 terms) and seeing which digits remain unchanged. If we required 20 d.p. of accuracy
then the solution is obtained once addition of successive terms causes no further change in the first 21
d.p. (21 to avoid rounding errors - it is prudent to obtain an extra digit of accuracy than required to avoid
rounding errors). This works because each successive term in the series is smaller than the proceeding
term - adding additional terms after the first few changes the value only slightly. The computed value
converges or settles down, to a value closer to pi with each term added.

Of course, if one wanted to beat the current world record by computing pi to trillions of decimal places,
then speed becomes a real issue. Such calculations are performed on supercomputers, and one would
probably use lower-level programming, avoiding the complexities (and niceties) of object-oriented
programming (OOP) which I have used, by using C++ or C, for extra speed (at the possible expense of
coding convenience). The last time I checked, the current world record was 10 trillion digits, calculated by
the Japanese mathematician Shigeru Kondo in 2011.

A very different approach to estimating pi, using probability theory, is illustrated in the article on Monte
Carlo methods.