Pi is a truly magical number that has a lot of people following it. I am not quite sure what is so fascinating with an irrational number that keeps repeating forever. From my point of view, I am interested in computing Pi, that is computing the value of Pi. Since Pi is an irrational number, it goes on forever. This means that any computation of Pi is just an approximation. If you compute 100 digits, I can compute 101 digits and be more accurate. Some people have set aside super computers to try to compute the most accurate Pi to date. Some extremes include Calculating 5 Trillion Digits of Pi. You can even find online a txt file containing 1 billion digits of Pi (Be warned, it will take a while to download the file and you can not open it with your regular notepad application). For me, it is interesting how you can compute Pi with a few simple line of Python.
You could always use the math.pi variable. It is included in the standard library and it should be used before you try to calculate it. In fact, we are going to use it to calculate our accuracy. To get started, let us take a look at a very straight forward approach of computing Pi. As usual, I will be using Python 2.7, same ideas might apply for different version. Most of the algorithms that we are going to use are direct implementations from the Pi WikiPedia page. Let us look at the following code:
import sys import math def main(argv): if len(argv) != 1: sys.exit('Usage: calc_pi.py <n>') print '\nComputing Pi v.01\n' a = 1.0 b = 1.0/math.sqrt(2) t = 1.0/4.0 p = 1.0 for i in range(int(sys.argv)): at = (a+b)/2 bt = math.sqrt(a*b) tt = t - p*(a-at)**2 pt = 2*p a = at;b = bt;t = tt;p = pt my_pi = (a+b)**2/(4*t) accuracy = 100*(math.pi-my_pi)/my_pi print "Pi is approximately: " + str(my_pi) print "Accuracy with math.pi: " + str(accuracy) if __name__ == "__main__": main(sys.argv[1:])
This is a very simple script, you can download, run, modify and share it as you like. You may see an output similler to the one bellow:
You will notice that one n is larger than 4 our approximated Pi and accuracy do not change anymore. We can only assume that the same will apply for larger values of n. So what are we to do? Lucky, there is more than one way to crack this egg. Using the Python Decimal Library we can compute Pi to a very high degree of accuracy. Let us see how the library works. The simple version is that it gives you a decimal number with more than 11 digits the normal Python float will give you. Here is an example of the Python Decimal Library:
Look at all those digits. Wait! We only entered 3.14, why are we getting all that junk? That is memory junk. In a nut shell, Python gave you the decimal number you wanted, plus a little extra. It won’t effect any calculation as long as the precision is less than when the junk numbers begin. You can specify how many digits you want exactly by setting the getcontext().prec to the number of digits you want. Let us try it out.
Sweet. Now let us try to use this to see if we can get a better approximation with our previous code. Now I normally am against using from library import *, but in this case it will make the code look nicer.
import sys import math from decimal import * def main(argv): if len(argv) != 1: sys.exit('Usage: calc_pi.py <n>') print '\nComputing Pi v.01\n' a = Decimal(1.0) b = Decimal(1.0/math.sqrt(2)) t = Decimal(1.0)/Decimal(4.0) p = Decimal(1.0) for i in range(int(sys.argv)): at = Decimal((a+b)/2) bt = Decimal(math.sqrt(a*b)) tt = Decimal(t - p*(a-at)**2) pt = Decimal(2*p) a = at;b = bt;t = tt;p = pt my_pi = (a+b)**2/(4*t) accuracy = 100*(Decimal(math.pi)-my_pi)/my_pi print "Pi is approximately: " + str(my_pi) print "Accuracy with math.pi: " + str(accuracy) if __name__ == "__main__": main(sys.argv[1:])
And the output:
Ok. We are more accurate, but it seems like there is some rounding off. We got the same accuracy from n=100 and from n=1,000. So now what? Well, now we come to ask for help from formulas. The way we calculated Pi until now is by summation of the parts. I found some code online from DAN on his article about Calculating Pi. He suggests we use the following 3 formulas:
Let us start with the Bailey–Borwein–Plouffe formula. It looks like this:
In code we can write it like this:
import sys import math from decimal import * def bbp(n): pi = Decimal(0) k = 0 while k < n: pi += (Decimal(1)/(16**k))*((Decimal(4)/(8*k+1))-(Decimal(2)/(8*k+4))-(Decimal(1)/(8*k+5))-(Decimal(1)/(8*k+6))) k += 1 return pi def main(argv): if len(argv) != 2: sys.exit('Usage: BaileyBorweinPlouffe.py <prec> <n>') getcontext().prec = (int(sys.argv)) my_pi = bbp(int(sys.argv)) accuracy = 100*(Decimal(math.pi)-my_pi)/my_pi print "Pi is approximately " + str(my_pi) print "Accuracy with math.pi: " + str(accuracy) if __name__ == "__main__": main(sys.argv[1:])
Set aside the “wrapping” code around it, the bbp(n) function is really what you want. The larger n you give it and the larger you set getcontext().prec the more accurate your calculation will be. Let us take a look at some exaction of the code:
That is plenty of digits. As you can tell, we are not more accurate than we were before. So we need to move on to our next formula, hoping to get a better accuracy, the Bellard’s formula. It looks like this:
We are going to only change our commutation formula, the rest of the code will remain the same. If you want you can download Bellard's formula in python. Let us take a look at bellards(n):
def bellard(n): pi = Decimal(0) k = 0 while k < n: pi += (Decimal(-1)**k/(1024**k))*( Decimal(256)/(10*k+1) + Decimal(1)/(10*k+9) - Decimal(64)/(10*k+3) - Decimal(32)/(4*k+1) - Decimal(4)/(10*k+5) - Decimal(4)/(10*k+7) -Decimal(1)/(4*k+3)) k += 1 pi = pi * 1/(2**6) return pi
Not , we get the same accuracy. Well, let us try the third formula, Chudnovsky algorithm, that looks like this:
Again, let us only look at the computation formula (assume we have a factorial formula). If you want you can download Chudnovsky formula in python.
Here is the code and output:
def chudnovsky(n): pi = Decimal(0) k = 0 while k < n: pi += (Decimal(-1)**k)*(Decimal(factorial(6*k))/((factorial(k)**3)*(factorial(3*k)))* (13591409+545140134*k)/(640320**(3*k))) k += 1 pi = pi * Decimal(10005).sqrt()/4270934400 pi = pi**(-1) return pi
So what is our conclusion? Fancy algorithms do not live up to the expectations in machine float worlds. I was really looking forward to a better accuracy than what we got when we used the summation formula. I guess that is too much to ask for. If you really need Pi, just use the math.pi variable. However, for fun and to test out how fast you computer really is, you can always try to calculate the first million digits of Pi, or maybe more...