Happy Pi Day!

Happy Pi Day Everyone! You must be wondering to yourself, “why is he celebrating pi day today?” Well, its because today is March 14, which is “3/14″, which is the first 3 digits of pi. I’ve been celebrating every pi day since 2006:

2006
2007
2008

And this year, I decided to talk about actually calculating pi. For thousands of years, people have been approximating  pi, which is nice, with some approximations getting pretty close, but these approximations don’t really give the most accurate values, which is why people started to calculate it.

I looked at a few algorithms, and I decided to use Gauss-Legendre’s algorithm, mainly because it was the simplest. I found that calculating pi on my regular machine to be quite difficult. I started off first in python:

def pi(a, b, t, p, count):
    if count > 0:
        a1 = (a + b)/2
        b1 = pow(a * b, 0.5)
        t1 = t - p*pow(a - a1, 2)
        p1 = 2 * p
        return pi(a1, b1, t1, p1, count - 1)
    else:
        return pow(a + b, 2)/(4*t)
 
a = 1.0
b = 1./pow(2, 0.5)
t = 1./4
p = 1.0
count = 100
print pi(a,b,t,p,count)

But I found that biggest value I could calculate was: 3.14159265359

I wanted something bigger. More precision! So I switched to Java, mainly because I would be able to use the “double” data type, which is a 64-bit floating point number. Unfortunately, I was some what limited as well:

import java.math.*;
 
public class pi {
    public static void main(String[] args) {
        double a = 1.0;
        double b = 1.0/Math.sqrt(2);
	double t = 1.0/4.0;
	double p = 1.0;
	int count = 1024;
 
	System.out.println(pi_calc(a, b, t, p, count));
    }
 
    private static double pi_calc(double a, double b, 
                                           double t, double p, int count){
	if (count > 0) {
	    double a1 = (a + b)/2;
	    double b1 = Math.sqrt(a * b);
	    double t1 = t - p*Math.pow(a - a1, 2);
	    double p1 = 2 * p;
 
            return pi_calc(a1, b1, t1, p1, count - 1);
        } else {
	    return Math.pow(a + b, 2)/(4*t);
        }
    }
}

The max count was only 1024, which would give me a value of 3.141592653589794. It is bigger then what python gave me, but not big enough! I mean, I wanted something that would give me 100000 digits of pi! So now I am looking into the java class, BigDecimal.

BigDecimal gives us a lot of useful methods, like add, subtract, divide, multiply, pow, etc… but its missing one crucial one, sqrt. Sure there’s a pow, but the pow method provided only accepts integers. So no “pow(0.5)” possible, which is why we need to write one.  While researching how to calculate sqrt, I stumbled upon this site, that had a sqrt calculating using BigDecimal! It was done through the use of Newton’s method:

public static BigDecimal sqrt(BigDecimal A, final int SCALE) {
    BigDecimal x0 = new BigDecimal("0");
    BigDecimal x1 = new BigDecimal(Math.sqrt(A.doubleValue()));
 
    while (!x0.equals(x1)) {
        x0 = x1;
        x1 = A.divide(x0, SCALE, BigDecimal.ROUND_HALF_UP);
        x1 = x1.add(x0);
        x1 = x1.divide(TWO, SCALE, BigDecimal.ROUND_HALF_UP);
    }
 
    return x1;
}

Now we have the tools to calculate pi! Coincidentally, that same page calculated pi as well, using a solution similar to mine. They used iteration instead of recursion, which is better in this scenario, since the max recursive depth wouldn’t be reached. So here’s our solutions combined:

import java.math.BigDecimal;
 
public class pi {
    final static BigDecimal TWO = new BigDecimal(2);
    final static BigDecimal FOUR = new BigDecimal(4);
 
    public static void main(String[] args) {
 
        int scale = 1024;
 
	BigDecimal a = new BigDecimal(1);
	BigDecimal b = a.divide(sqrt(TWO, scale), scale, BigDecimal.ROUND_HALF_UP);
	BigDecimal t = new BigDecimal(0.25);
	BigDecimal p = a;
	BigDecimal y;
 
	while (!a.equals(b)) {
            y = a;
            a = a.add(b).divide(TWO, 
                                      scale, BigDecimal.ROUND_HALF_UP);
            b = sqrt(b.multiply(y), scale);
            t = t.subtract(p.multiply(y.subtract(a).multiply(y.subtract(a))));
            p = p.multiply(TWO);
        }
 
	BigDecimal result = a.add(b).multiply(a.add(b)).divide(t.multiply(new BigDecimal(4)),
	                                                scale, BigDecimal.ROUND_HALF_UP);
 
        System.out.println(result);
	System.out.println(result.precision());
    }
}

So what’s the precision? Well see that scale variable? Yes, that’s how many digits I have:

3.14159265358979323846264338327950288419716939937510582097494459230
7816406286208998628034825342117067982148086513282306647093844609550
5822317253594081284811174502841027019385211055596446229489549303819
6442881097566593344612847564823378678316527120190914564856692346034
8610454326648213393607260249141273724587006606315588174881520920962
8292540917153643678925903600113305305488204665213841469519415116094
3305727036575959195309218611738193261179310511854807446237996274956
7351885752724891227938183011949129833673362440656643086021394946395
2247371907021798609437027705392171762931767523846748184676694051320
0056812714526356082778577134275778960917363717872146844090122495343
0146549585371050792279689258923542019956112129021960864034418159813
6297747713099605187072113499999983729780499510597317328160963185950
2445945534690830264252230825334468503526193118817101000313783875288
6587533208381420617177669147303598253490428755468731159562863882353
7875937519577818577805321712268066130019278766111959092164201989380
952572010654858632805

Precision: 1025

Enjoy your pi!