%pylab inline
import scipy
Populating the interactive namespace from numpy and matplotlib
def monte_carlo_pi(N):
"""
Percent of area of square from (0,0) to (1,1)
that is inside x^2 + y^2 < 1 is pi/4.
"""
x = scipy.random.rand(N)
y = scipy.random.rand(N)
pi = 4.0*sum(x*x+y*y < 1)/N
pi_err = abs(pi - scipy.pi)
return pi, pi_err
def leibniz_pi(N):
"""
pi/4 = 1 - 1/3 + 1/5 - 1/7 + 1/9 - ...
Obtained by Taylor expansion of arctan (arctan(x) = x - x^3/3 + x^5/5 - x^7/7 + x^9/9 + ... )
Also, noting arctan(1) = pi/4
"""
pi = 4*sum((-1)**i/(2*i+1.0) for i in range(N))
pi_err = abs(pi - scipy.pi)
return pi, pi_err
def machin_pi(N):
"""
pi/4 = 4 arctan(1/5) - arctan(1/239)
"""
def arctan(x, N):
return sum((-1)**i *(x**(2*i+1.0))/(2*i+1.0) for i in range(N))
pi = 16*arctan(1/5., N//2) - 4 * arctan(1/239.0, N//2)
pi_err = abs(pi - scipy.pi)
return pi, pi_err
z = scipy.arange(6,100)
z_exp = scipy.floor((1.2)**z)
data = {}
for f in (monte_carlo_pi, leibniz_pi, machin_pi):
data[f] = []
for z in z_exp:
data[f].append(f(int(z)))
mc, mc_err = zip(*data[monte_carlo_pi])
leibniz, leibniz_err = zip(*data[leibniz_pi])
machin, machin_err = zip(*data[machin_pi])
plot(z_exp,mc, 'r', label="Monte Carlo")
plot(z_exp, leibniz, 'b', label="Leibniz")
plot(z_exp, machin, 'g', label="Machin")
xscale('log')
axis([1,10**7,0,5])
legend(loc='lower right')
<matplotlib.legend.Legend at 0x107d5d1d0>
plot(z_exp,mc_err, 'r', label = "MC Error")
plot(z_exp, leibniz_err, 'b', label="Leibniz Error")
plot(z_exp, machin_err, 'g', label="Machin Error")
xscale('log')
yscale('log')
legend(loc='lower right')
<matplotlib.legend.Legend at 0x108f86fd0>
import timeit
print "Monte Carlo Time:", timeit.repeat("monte_carlo_pi(100000)", setup="""def monte_carlo_pi(N):
import scipy
x = scipy.random.rand(N)
y = scipy.random.rand(N)
pi = 4.0*sum(x*x+y*y < 1)/N
pi_err = abs(pi - scipy.pi)
return pi, pi_err""", number=10, repeat=3)
print "Leibniz Time", timeit.repeat("leibniz_pi(100000)", setup="""def leibniz_pi(N):
pi = 4*sum((-1)**i/(2*i+1.0) for i in range(N))
pi_err = abs(pi - 3.141592653589793238)
return pi, pi_err""", number=10, repeat=3)
print "Leibniz (no -1^i by expanding out two terms) Time", timeit.repeat("leibniz_pi(100000)", setup="""def leibniz_pi(N):
pi = 4*sum(1.0/(4*i+1.0) - 1.0/(4*i+3.0) for i in range(N//2))
pi_err = abs(pi - 3.141592653589793238)
return pi, pi_err""", number=10, repeat=3)
print "Machin Time", timeit.repeat("machin_pi(100000)", setup="""def machin_pi(N):
def arctan(x, N):
return sum((-1)**i *(x**(2*i+1.0))/(2*i+1.0) for i in range(N))
pi = 16*arctan(1/5., N//2) - 4 * arctan(1/239.0, N//2)
pi_err = abs(pi - 3.141592653589793238)
return pi, pi_err""", number=10, repeat=3)
print "Machin (N=22) Time", timeit.repeat("machin_pi(22)", setup="""def machin_pi(N):
def arctan(x, N):
return sum((-1)**i *(x**(2*i+1.0))/(2*i+1.0) for i in range(N))
pi = 16*arctan(1/5., N//2) - 4 * arctan(1/239.0, N//2)
pi_err = abs(pi - 3.141592653589793238)
return pi, pi_err""", number=10, repeat=3)
Monte Carlo Time: [2.4383060932159424, 2.4071178436279297, 2.343600034713745] Leibniz Time [0.44296717643737793, 0.43317699432373047, 0.45170092582702637] Leibniz (no -1^i by expanding out two terms) Time [0.12793684005737305, 0.1200859546661377, 0.12012195587158203] Machin Time [0.5834639072418213, 0.5719470977783203, 0.5927698612213135] Machin (N=22) Time [0.00010800361633300781, 0.00011205673217773438, 0.00013113021850585938]
Note, MC is slowest and least accurate method.
Leibniz and Machin take comparable times until you recognize that you don't need to calculate $(-1)^i$ if you unroll two terms in the summation (one positive, one negative), at which point Leibniz is more than 4 times faster (by skipping exponentiation step) on a term by term basis.
However, Machin is ultimately more than 1000 times faster since you don't need more than 22 terms (11 in each arctan) when you use 64-bit floating point numbers to get to the limit of floating point accuracy. It also gets to an accuracy of about 10^-15 while the 100000 term Leibniz series only gets to an accuracy of only 10^-5, so the Machin method at 22 terms is also 10 billion times more accurate in addition to being more than 1000 times faster.
print "Machin(22):", machin_pi(22)
print "Leibniz(100000):", leibniz_pi(100000)
Machin(22): (3.141592653589794, 8.881784197001252e-16) Leibniz(100000): (3.1415826535897198, 1.0000000073340232e-05)
Note, you could use an arbitrary precision floating point package like the following to increase the accuracy of the Machin method:
from mpmath import mp
def machin_mp_pi(N):
"""
pi/4 = 4 arctan(1/5) - arctan(1/239)
"""
mp.dps = N if N > 25 else 25
def arctan(x, N):
return sum((x**(4*i+1.0))/(4*i+1.0) - (x**(4*i+3.0))/(4*i+3.0) for i in range(N//2))
pi = 16*arctan(mp.mpf("1/5"), N//2) - 4 * arctan(mp.mpf("1/239"), N//2)
pi_err = float(abs(pi - mp.pi))
return pi, pi_err
for d in [10,20,30,40,50,60,70,80,90,100,200,300,400]:
print "N: %d, Error: %g" % (d, machin_mp_pi(d)[1])
N: 10, Error: 8.81408e-07 N: 20, Error: 1.54155e-15 N: 30, Error: 2.85522e-21 N: 40, Error: 8.26631e-30 N: 50, Error: 1.77018e-35 N: 60, Error: 5.82261e-44 N: 70, Error: 1.31758e-49 N: 80, Error: 4.59657e-58 N: 90, Error: 1.07086e-63 N: 100, Error: 3.86473e-72 N: 200, Error: 2.46085e-142 N: 300, Error: 2.08286e-212 N: 400, Error: 1.98177e-282