Saturday, October 11, 2008

Calculation of Pi

This time I tried to answer a question - what is the best method of calculating Pi to a required accuracy - for example,- 100 or 1000 or 10 000 decimal digits ?
So I chose 4 Pi formulas for testing:
1. Leibniz
2. Madhava transformed series
3. Machin like formula
4. Gauss algorithm
The idea for experiment - is to find out which algorithm gives more decimal digits of Pi in same number of iterations/terms. So measurement algorithm is this:
1. Calculate correct digits of Pi given from every Pi formula in each iteration.
2. Plot diagram - correct Pi decimal digits amount VS iteration number
(For this to work - you must have Matplotlib module).
Python code doing this stuff:



import inspect, time, sys
from decimal import Decimal, getcontext
from pylab import plot, show, xlabel, ylabel, title, grid, legend

def correct_digits(calcpi):
pi = "3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989"
gd = 1000
for i,s in enumerate(zip(pi,str(calcpi))):
if s[0] != s[1]:
gd = i - 2
if gd < 0: gd = 0
break
return gd

def show_results(res):
x=range(len(res[0][2]))
plpar = ""
legpar = ""
for i in range(len(res)):
if plpar: plpar +=","
if legpar: legpar +=","
plpar += "x, res["+str(i)+"][2]"
legpar += "res["+str(i)+"][0][3:] +' ('+ str(int(res["+str(i)+"][1])) + 'ms)'"
plpar = "plot("+plpar+", linewidth=2.0)"
legpar = "legend(("+legpar+"), loc=0)"
exec(plpar)
exec(legpar)
xlabel('Iterations', fontsize = 15)
ylabel('Correct digits', fontsize = 15)
title("Convergence of some $\pi$ formulas",fontsize=20)
grid()
show()

def test_pi_formulas(pi_func_list):
res = []
terms = 6
print 'Calculating %d iterations...' % terms
getcontext().prec = 1002
for pif in pi_func_list:
t1 = time.time()*1000
name,r = pif(terms)
t2 = time.time()*1000
res.append((name,t2-t1,r))
show_results(res)

def pi_Leibniz(terms):
"""
http://en.wikipedia.org/wiki/Leibniz_formula_for_pi
"""
st = inspect.stack()[0][3]
ra = Decimal(0)
ans = []

for n in range(terms):
ra += Decimal(4*(-1)**n)/Decimal(2*n+1)
ans.append(correct_digits(ra))

return st,ans

def pi_Madhava(terms):
"""
http://en.wikipedia.org/wiki/Madhava_of_Sangamagrama#The_value_of_.CF.80_.28pi.29
"""
st = inspect.stack()[0][3]
ra = Decimal(0)
k = Decimal(12)**Decimal('0.5')
ans = []

for n in range(terms):
ra += Decimal((-1)**n)/Decimal((2*n+1)*(3**n))
ans.append(correct_digits(k*ra))

return st, ans

def pi_Machin(terms):
"""
http://en.wikipedia.org/wiki/Machin-like_formula
"""
st = inspect.stack()[0][3]
c1 = Decimal(1)/Decimal(5)
c2 = Decimal(1)/Decimal(239)
a = Decimal(4)
atan1 = Decimal(0)
atan2 = Decimal(0)
ra = Decimal(0)
ans = []

for n in range(terms):
atan1 += Decimal((-1)**n) * (c1**Decimal(2*n+1))/Decimal(2*n+1)
atan2 += Decimal((-1)**n) * (c2**Decimal(2*n+1))/Decimal(2*n+1)
ra = a * (a*atan1 - atan2)
ans.append(correct_digits(ra))

return st, ans

def pi_Gauss(terms):
"""
http://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_algorithm
"""
st = inspect.stack()[0][3]
one = Decimal(1)
ra = Decimal(0)
val = (one,one/(Decimal(2)**Decimal('0.5')),one/Decimal(4),one)
ans = []

for n in range(terms):
an, bn, tn, pn = val
an1 = (an+bn)/Decimal(2)
bn1 = (an*bn)**Decimal('0.5')
tn1 = tn - pn * (an - an1)**Decimal(2)
pn1 = Decimal(2) * pn
val = (an1, bn1, tn1, pn1)
ra = ((an+bn)**Decimal(2))/(Decimal(4)*tn)
ans.append(correct_digits(ra))

return st, ans

if __name__ == "__main__":
mod = sys.modules[__name__]
pifunc = []
for key, value in mod.__dict__.items():
if str(value)[:13] == '<function pi_':
pifunc.append(value)
test_pi_formulas(pifunc)



BTW this code is very scalable - if there is a need to test additional Pi formula - you only need to add new Pi function code like this:
"def pi_new_required_function(terms):
[...]". (Function name MUST begin with "pi_"). The rest program will do it for you - measure execution times, measure decimal digits correctness, plot results...
This is thanks to Python reflection used in program -> program scans itself for ALL pi functions and dynamically invokes them.

So after execution of this test we get such graph:


Conclusion
From graph can be seen that fastest converges Gauss Pi algorithm (from the chosen set of Pi formulas). So if you need for some reason to calculate Pi in desired accuracy (desired decimal places amount) I suggest Gauss algorithm. Because you will get result in less iterations / CPU cycles compared with other algorithms from the set.

Have fun!

No comments:

Post a Comment

Comment will be posted after comment moderation.
Thank you for your appreciation.