Mathematics

Computing Pi with Python

May 28, 2013

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.

wpid-200px-Innovation_and_Technology_Commission.svg-2013-05-28-12-54.png

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[1])):
		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:

wpid-simple_python_pi_calculation-2013-05-28-12-54.png
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:

wpid-python_decimal_example-2013-05-28-12-54.png

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.

wpid-python_decimal_example_6_digits-2013-05-28-12-54.png

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[1])):
		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:

wpid-python_pi_accuracy-2013-05-28-12-54.png

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:

wpid-48f7653d58f4ad747327d271ed789415-2013-05-28-12-54.png

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 &lt; 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[1]))
	my_pi = bbp(int(sys.argv[2]))
	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:

wpid-python_BaileyBorweinPlouffe-2013-05-28-12-54.png

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:

wpid-dbf2d4355c108f6b3388985be4976799-2013-05-28-12-54.png

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 &lt; 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

The output:

wpid-bellards_pi_run-2013-05-28-12-54.png

Not , we get the same accuracy. Well, let us try the third formula, Chudnovsky algorithm, that looks like this:

wpid-826dc7788dba249ee86fc0135e06b035-2013-05-28-12-54.png

 

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 &lt; 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

wpid-chudnovsky_pyhton-2013-05-28-12-54.png

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

If you enjoyed this post, please consider leaving a comment or subscribing to the RSS feed to have future articles delivered to your feed reader.

You Might Also Like

  • Person

    Hi, I believe the reason Decimal(3.14) gives you "junk" is that the
    input is first interpreted as a float (which won't be 3.14 exactly), not
    memory junk. Decimal(0.125) is OK, presumably because it is exactly expressible in binary floating point. Decimal('3.14') is also OK since it circumvents that issue.

    • CptDeadBones

      That is precisely why. I just did not want to get into that level of technicality. It has to do with machine representations of numbers. Maybe I should address that in a future post. Thanks for pointing that out.
      Captain DeadBones
      thelivingpearl.com

  • ibex

    Nice job with wrong conclusion. The math.pi only has 16 digital. So the last comparison with it as a standard is meanless...

    • CptDeadBones

      Thanks for the comment!

      You are right about Math.pi only being 16 digits long. However, since we are comparing all the results to it it does not. Notice that I did not say that the error is off of true value of Pi, rather Math.pi value. Since all results compare to it and all accuracy results are the same, it is only logical that the conclusion stands.
      Captain DeadBones
      thelivingpearl.com

  • obitus

    Great post, very well documented. I like how you show the contrast in different algorithms.

    • CptDeadBones

      Thanks!

      Captain DeadBones
      thelivingpearl.com

  • Mohit S Gujarathi

    Hi, i am beginner in python .
    Could you please explain some basics from your code in 1 line each ?
    Would be appreciated a lot.

    1)what does this mean
    "" if len(argv) != 1:
    sys.exit('Usage: calc_pi.py ') ""
    in the very first program under the main function?

    2) Why cant we write """ from import decimal * """ as simply " import decimal " ?

    PS: I know very well my questions are lame but I am teaching myself coding without any collg education and only using the internet. So your help would be greatly appreciated.

    • CptDeadBones

      Hi Mohit!

      I am glad you asked these questions, there are always people who do not understand, but are too afraid to ask. Please feel free to ask any question you have that relates to my posts or Python and I will do my best to answer.

      1) The first if statement checks to see how many argument got passed in when the program was executed. Normally when we run a python script we will type:

      python example.py

      This program will have 0 command line arguments. If we wanted to pass in some variables, we can pass in an argument after the file name, like so:

      python example.py 5

      In this case the number 5 will be passed in. Going back to our if statement. In my example we want the user to pass in a number when excuating the program. If a number is not provided the program terminates. That's all.

      2) Either import statement works. The reason I chose that one is becuase otherwise I would have to write "decimal.Decimal()" instead of just "Decimal()" and that will make line 9 a very long statment.

      I hope that clears the air a little. Let me know if you have more questions!

      • Mohit S Gujarathi

        Thanks a lot . It became quite clear now.

        But I do require a lot of your help.
        I have finished this course from MIT OCW in a span of a week-http://ocw.mit.edu/courses/electrical-engineering-and-computer-science/6-00-introduction-to-computer-science-and-programming-fall-2008/video-lectures/

        Considering that as level 1, my aim is now to get to get level 2 and level 3.
        But I dont know what exactly constitutes those two levels.
        Which projects should I start working on to gain experience ?

        What I need is some definite roadmap on how to explore and gain skills. I would just need pointers and the rest I ll google.

        • CptDeadBones

          I am going to have to take a closer look at what exactly level 1, 2 and 3 are. Can you please send me an email so we can discuss this in private?

      • Mohit S Gujarathi

        Also was there a special reason behind writing the last if statement ==>

        if __name__ == "__main__":

        main(sys.argv[1:])

        Why didnt you simply call the function like this ==>

        main(sys.argv[1:])

        • CptDeadBones

          Nothing special. Again, both versions will work just fine.

  • Pingback: Finding yourself in pi | Dead End Math()

  • Pingback: Google()

  • rilian

    I get a syntax error w/ the bbp function in line: "while k < n:" when simply copying and pasting your code and running it. After looking at it a bit and doing a quick google search, it dawned on me that < is a replacement for '<'. Is this webpage replacing the character errantly? When I replace '<' with '<' in the code, everything runs as you describe...

    This also appears in the code for chudnovsky in my browser. I use chrome.

  • u8y7541

    Well... I used the algebraic circle formula to calculate all the points on a circle, then calculate the distance between each one and its neighbor and add everything up to get an approx. of the circumference. Then I divide by diameter. I got pi up to 11 digits...

  • u8y7541

    Nice, but this kind of seems completely useless. If you're approximating the number to "math.pi," why can't you just use math.pi in the first place? So you're not really computing pi. Real "computing pi" would be if you computed pi from SCRATCH.