## 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
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, sysfrom decimal import Decimal, getcontextfrom pylab import plot,  show,  xlabel,  ylabel,  title,  grid,  legenddef 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 gddef 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,ansdef 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, ansdef 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, ansdef 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, ansif __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!