CODEWITHSUNDEEP
pythonpython-3.xtrigonometry
Johnathan
Johnathan

Reputation: 1907

Python: estimate Pi with trig functions as efficiently as possible

I have an assignment where I need to approximate Pi in a computationally efficient manner. Here is my strategy: I use a unit circle, the angle bisector of an isoceles triangle, and the definition of sin. I drew a diagram:

enter image description here

For example, if I want to use an hexagon (6 points/6 sides), I simply need to compute a:(0.5*sin(2*pi/2*x) and multiply it by (2*x). Finally, since Pi = Circumference/Diameter, then my approximation of Pi = polygon perimeter (since Diameter = 1).

Essentially:

from math import sin, pi
def computePi(x):    #x: number of points desired
    p = x*sin(pi/x)
    print(p)

computePi(10000)
3.141592601912665

It works, and I think it's as efficient as it gets, no? Thank you for your time!

EDIT: to avoid circularity, I redid it using Archimedes algorithm using only the Pythagorean theroem:

enter image description here

Code:

from math import sqrt

def approxPi(x):                  #x: number of times you want to recursively apply Archmidedes' algorithm
    s = 1                         #Unit circle
    a = None; b = None;   
    for i in range(x):
        a = sqrt(1 - (s/2)**2)
        b = 1 - a
        print('The approximate value of Pi using a {:5g}-sided polygon is {:1.8f}'.format(6*2**(i),(s*6*2**(i))/2))
        s = sqrt(b**2 + (s/2)**2)

Upvotes: 5

Views: 2052

Answers (3)

John Coleman
John Coleman

Reputation: 51998

A fun albeit not very efficient solution is to use Euler's solution of the Basel Problem:

from math import sqrt

def psum(n):
    return sum(1/k**2 for k in range(1,n+1))

def approxPi(n):
    s = psum(n)
    return sqrt(6*s)

For example, 

>>> approxPi(100000)
3.141583104326456

As I said, not very efficient. On the other hand, there is clearly no subtle circularity. Many other series are known to either converge to pi or to converge to a value from which pi can be easily computed, and many of these other series converge much more rapidly.

On edit: @Simon 's suggestion of using the Gauss-Legendre algorithm, together with the module decimal, leads to the following code (which returns the result as a string):

import decimal
from decimal import Decimal as d

def approxPi(n):
    eps = 1/d(10**n)
    decimal.getcontext().prec = 3*n #probably overkill, but need room for products
    a = d(1)
    b = 1/d(2).sqrt()
    t = 1/d(4)
    p = d(1)
    dif = a-b
    if dif < 0: dif = -dif
    i = 1
    while dif >= eps:
        a1 = (a+b)/2
        b1 = a*b
        b1 = b1.sqrt()
        t1 = t - p*(a - a1)**2
        p1 = 2*p
        a,b,t,p = a1,b1,t1,p1
        dif = a1-b1
        if dif < 0: dif = -dif
    pi = (a + b)**2/(4*t)
    return str(pi)[:n+2]

For example,

>>> approxPi(1000)
'3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989'

With agrees with this.

The above took less than a second. It took several seconds for 10,000. It would be interesting to see how long it would take to get 1,000,000 digits in Python with this.

Upvotes: 2

JonahR
JonahR

Reputation: 104

Here is the code for your problem:

from math import radians, sin


def computePi(n):
    p = n * (sin(radians((360/(2*n)))))
    print(p)
computePi(1000)

The theory behind this code is explained in this thread: https://math.stackexchange.com/questions/588141/how-is-the-value-of-pi-pi-actually-calculated

Upvotes: 1

Rory Daulton
Rory Daulton

Reputation: 22544

Even better is

print(4 * math.atan(1))

This does not use pi in any obvious way in the calculation (though as @Jean-FrançoisFabre comments, pi is probably used in the function definition), and in addition to the trig function it has just one simple multiplication. Of course, there is also

print(2 * math.acos(0))

and

print(2 * math.asin(1))

Upvotes: 5

Related Questions