The value of π from Python operations in the implementation of scientific computation in Python

Source: Internet
Author: User
Tags value of pi
π is a truly magical number followed by countless people. I am not very clear about the fascinating point of an irrational number that is always repeated. In my opinion, I am happy to calculate π, which is the value of pi. Because π is an irrational number, it is infinite. This means that any calculation of π is just an approximation. If you calculate 100 bits, I can calculate 101 bits and be more precise. So far, some people have chosen supercomputers to try to calculate the most accurate π. Some extrema include the 500 million bits that calculate π. You can even find a 10 billion-bit text file containing π from the Internet (note!). It may take a while to download this file, and it won't open with the Notepad app you normally use. )。 For me, how to calculate π with a few lines of Python is my interest.
You can always use the Math.PI variable. It is included in the standard library, and you should use it before you try to calculate it yourself. In fact, we will use it to calculate the accuracy. As a starting point, let's look at a very straightforward method of calculating pi. As usual, I will use Python 2.7, and the same ideas and code might apply to different versions. Most of the algorithms we will be using are from the Pi WikiPedia page and implemented. Let's take a look at the following code:

Importsysimportmath Defmain (argv):   iflen (argv)!=1:    sys.exit (' Usage:calc_pi.py 
 
  
   
  ')   print ' \ncomputing Pi v.01\n '     a=1.0  b=1.0/math.sqrt (2)  t=1.0/4.0  p=1.0       foriinrange (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 that you can download, run, modify, and share freely with others. You can see output similar to the following:

You will find that even though N is greater than 4, we are approaching Pi accuracy without much improvement. We can guess that even if the value of n is large, the same thing (pi approximation accuracy does not improve) will still occur. Fortunately, there is more than one way to uncover this mystery. Using the Python decimal (decimal) library, we can get a higher precision value to approximate the pi. Let's look at how the library function is used. With this simplified version, you can get more than 11 digits of the number typically given by the small Python floating point precision. Here is an example from the Python Decimal library:

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

See these numbers. Wrong! We entered only 3.14, why did we get some rubbish (junk)? This is memory garbage (junk). Simply put, Python gives you the decimal number you want, plus a little extra value. As long as the precision is less than the garbage count, it does not affect any calculations. By setting GetContext (). Prec you can get to the number of digits you want. Let's try.

See these numbers. Wrong! We entered only 3.14, why did we get some rubbish (junk)? This is memory garbage (junk). Simply put, Python gives you the decimal number you want, plus a little extra value. As long as the precision is less than the garbage count, it does not affect any calculations. By setting GetContext (). Prec you can get to the number of digits you want. Let's try.

Very good. Now let's try to use this to see if we can get a better approximation of our previous code. Now, I usually object to using "From library import *", but in this case it will make the code look more beautiful.

importsysimportmathfromdecimalimport* Defmain (argv):   iflen (argv)!=1:    sys.exit (' Usage:calc_pi.py 
 
  
   
  ')   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)       foriinrange (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:])
 
  


Output Result:

All right. We're more accurate, but it looks like there's some rounding. from n = 100 and n = 1000, we have the same precision. Now what? OK, now let's turn to the formula. So far, the way we calculate pi is by adding a few parts together. I found some code from Dan's article on calculating Pi. He suggested that we use the following 3 formulas:

Bailey–borwein–plouffe formula
The formula of Bellard
Chudnovsky algorithm

Let's start with the Bailey–borwein–plouffe formula. It looks like this:

In the code we can write it this way:

Import sysimport mathfrom 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 
 
   
  
   
    
   ')       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:])
  
   
 
  


Throw away the "wrapper" code, the BBP (N) function is what you really want. You give it the greater the N and give the GetContext (). Prec the larger the value, the more accurate the calculation will be. Let's take a look at some of the code results:

This has many digit digits. As you can see, we are not more accurate than before. So we need to go forward to the next formula, Beira the formula, hoping to get better accuracy. It looks like this:

We will only change our transformation formula, and the rest of the code will remain unchanged. Click here to download the Beira formula implemented by Python. Let's 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 (+)/(10*k+3)-decimal (+)/(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

Oh, no, we got the same precision. Well, let's try the third formula, the Chudnovsky algorithm, which looks like this:

Again, let's take a look at this formula (suppose we have a factorial formula). Click here to download the Chudnovsky formula implemented in Python.

Here are the program and output results:

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 do we have to conclude? Fancy algorithms do not make the machine floating-point world reach a higher standard. I really look forward to having a better precision than we can get with the summation formula. I guess that's too much to ask. If you really need pi, just use the Math.PI variable. However, as fun and test how fast your computer can really be, you can always try the first one to calculate the PI million dollar or more bits that are several.

  • Related Article

    Contact Us

    The content source of this page is from Internet, which doesn't represent Alibaba Cloud's opinion; products and services mentioned on that page don't have any relationship with Alibaba Cloud. If the content of the page makes you feel confusing, please write us an email, we will handle the problem within 5 days after receiving your email.

    If you find any instances of plagiarism from the community, please send an email to: info-contact@alibabacloud.com and provide relevant evidence. A staff member will contact you within 5 working days.

    A Free Trial That Lets You Build Big!

    Start building with 50+ products and up to 12 months usage for Elastic Compute Service

    • Sales Support

      1 on 1 presale consultation

    • After-Sales Support

      24/7 Technical Support 6 Free Tickets per Quarter Faster Response

    • Alibaba Cloud offers highly flexible support services tailored to meet your exact needs.