Table of Contents¶
Python Pi Day 2017¶
This is heavily inspired by what I did two years ago, see this page cs101/hackhaton/14_03/2015 on my website.
Today is Pi Day 2017, the day celebrating the number $\pi$. For more details on this number, see this Wikipedia page.
Let us use this occasion to showcase a few different approaches to compute the digits of the number $\pi$. I will use the Python programming language, and you are reading a Jupyter notebook.
Computing a lot of digits of $\pi$?¶
- What to do ? I will present and implement some methods that can be used to compute the first digits of the irrational number $\pi$.
- How ? One method is based on random numbers, but all the other are simple (or not so simple) iterative algorithm: the more steps you compute, the more digits you will have!
- How to compare / assess the result ? It is simple: the more digits you got, the better. We will also test the different methods implemented, and for the most efficient, see how fast it is to go up-to $140000$ digits.
The simple goal is to write a small function that computes digits of pi, as fast as possible, and find the 10 digits from position 140317 to 140327! (that's the challenge I posted on Facebook)
Three approximations, using fractions, were known from a very long time (Aristote used $\frac{355}{113}$ !). The first three approximations of pi are:
print(" pi ~= 3.14 (two first digits).")
print(" pi ~= 22/7 = {} (two first digits).".format(22.0 / 7.0))
print(" pi ~= 355/113 = {} (six first digits).".format(355.0 / 113.0))
pi ~= 3.14 (two first digits). pi ~= 22/7 = 3.142857142857143 (two first digits). pi ~= 355/113 = 3.1415929203539825 (six first digits).
This method is extremely limited, and will not give you more than 13 correct digits, as math.pi is stored as a float number (limited precision).
def mathpi():
from math import pi
return pi
print("First method: using math.pi gives pi ~= {:.17f} (17 digits are displayed here).".format(mathpi()))
First method: using math.pi gives pi ~= 3.14159265358979312 (17 digits are displayed here).
If we know the digits, we can directly print them:
from decimal import Decimal
bigpi = Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679')
print("The first 100 digits of pi are {}.".format(bigpi))
The first 100 digits of pi are 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679.
A simple Monte-Carlo method¶
A simple Monte Carlo method for computing $\pi$ is to draw a circle inscribed in a square, and randomly place dots in the square. The ratio of dots inside the circle to the total number of dots will approximately equal $\pi / 4$, which is also the ratio of the area of the circle by the area of the square:
In Python, such simulation can basically be implemented like this:
from random import uniform
def montecarlo_pi(nbPoints=10000):
"""Returns a probabilist estimate of pi, as a float number."""
nbInside = 0
# we pick a certain number of points (nbPoints)
for i in range(nbPoints):
x = uniform(0, 1)
y = uniform(0, 1)
# (x, y) is now a random point in the square [0, 1] × [0, 1]
if (x**2 + y**2) < 1:
# This point (x, y) is inside the circle C(0, 1)
nbInside += 1
return 4 * float(nbInside) / floor(nbPoints)
print("The simple Monte-Carlo method with {} random points gave pi = {}".format(10000, montecarlo_pi(10000)))
The simple Monte-Carlo method with 10000 random points gave pi = 3.138
It is an interesting method, but it is just too limited for approximating digits of $\pi$.
The previous two methods are quite limited, and not efficient enough if you want more than 13 digits (resp. 4 digits for the Monte-Carlo method).
$100$ first digits of $\pi$¶
$\pi \simeq 3.1415926535 ~ 8979323846 ~ 2643383279 ~ 5028841971 \\\\ 6939937510 ~ 5820974944 ~ 5923078164 ~ 9862803482 ~ 53421170679$ when computed to the first $100$ digits.
Can you compute up to $1000$ digits? Up to $10000$ digits? Up to $100000$ digits? Up to 1 million digits?
Using mpmath¶
This is a simple and lazy method, using the mpmath module.
MPmath is a Python library for arbitrary-precision floating-point arithmetic (Multi-Precision), and it has a builtin highly-optimized algorithm to compute digits of $\pi$.
It works really fine up-to 1000000 digits (56 ms), from 1 million digits to be printed, printing them starts to get too time consuming (the IDE or the system might freeze).
import mpmath
# from sympy import mpmath # on older sympy versions
mp = mpmath.mp
We can arbitrarily set the precision, with the constant mp.dps (digit numbers).
mp.dps = 1000 # number of digits
my_pi = mp.pi # Gives pi to a thousand places
print("A lazy method using the mpmath module:\npi is approximatly {} (with {} digits).".format(my_pi, mp.dps))
A lazy method using the mpmath module: pi is approximatly 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420198 (with 1000 digits).
Let save it for further comparison of simpler methods.
mp.dps = 100000 # number of digits
len(str(mp.pi))
mpmath_pi = Decimal(str(mp.pi))
100001
We can solve the initial challenge easily:
mp.dps = 140330
print(str(mp.pi)[2:][140317:140317+10])
9341076406
And it will most probably be the quickest method presented here:
%timeit mp.dps=140330;print(str(mp.pi)[2:][140317:140317+10])
9341076406 9341076406 9341076406 9341076406 1 loop, best of 3: 236 ms per loop
Asking for $10$ times more digits take about $100$ more of time (that's a bad news).
%timeit mp.dps=1403230;print(str(mp.pi)[2:][1403217:1403217+10])
7966916520 7966916520 7966916520 7966916520 The slowest run took 6.34 times longer than the fastest. This could mean that an intermediate result is being cached. 1 loop, best of 3: 20.4 s per loop
The Gauss–Legendre iterative algorithm¶
More details can be found on this page.
The first method given here is explained in detail. This algorithm is called the Gauss-Legendre method, and for example it was used to compute the first $206 158 430 000$ decimal digits of $\pi$ on September 18th to 20th, $1999$.
It is a very simply iterative algorithm (ie. step by step, you update the variables, as long as you want):
First, start with 4 variables $a_0, b_0, t_0, p_0$, and their initial values are $a_0 = 1, b_0 = 1/\sqrt{2}, t_0 = 1/4, p_0 = 1$.
Then, update the variables (or create 4 new ones $a_{n+1}, b_{n+1}, t_{n+1}, p_{n+1}$) by repeating the following instructions (in this order) until the difference of $a_{n}$ and $b_{n}$, is within the desired accuracy:
- $a_{n+1} = \frac{a_n + b_n}{2}$,
- $b_{n+1} = \sqrt{a_n \times b_n}$,
- $t_{n+1} = t_n - p_n (a_n - a_{n+1})^2$,
- $p_{n+1} = 2 p_n$.
Finally, the desired approximation of $\pi$ is: $$\pi \simeq \frac{(a_{n+1} + b_{n+1})^2}{4 t_{n+1}}.$$
Using float numbers¶
The first three iterations give (approximations given up to and including the first incorrect digit):
3.140 …
3.14159264 …
3.1415926535897932382 …
The algorithm has second-order convergent nature, which essentially means that the number of correct digits doubles with each step of the algorithm. Try to see how far it can go in less than 120 seconds (print the approximation of $\pi$ after every step, and stop/kill the program after 2 minutes).
This method is based on computing the limits of the arithmetic–geometric mean of some well-chosen number (see on Wikipédia for more mathematical details).
import math
def gauss_legendre_1(max_step):
"""Float number implementation of the Gauss-Legendre algorithm, for max_step steps."""
a = 1.
b = 1./math.sqrt(2)
t = 1./4.0
p = 1.
for i in range(max_step):
at = (a + b) / 2.0
bt = math.sqrt(a*b)
tt = t - p*(a-at)**2
pt = 2.0 * p
a, b, t, p = at, bt, tt, pt
my_pi = ((a+b)**2)/(4.0*t)
return my_pi
my_pi = gauss_legendre_1(4)
my_pi
print("pi is approximately: {:.15f} (as a float number, precision is limited).".format(my_pi))
accuracy = 100*abs(math.pi - my_pi)/math.pi
print("Accuracy % with math.pi: {:.4g}".format(accuracy))
accuracy = 100*abs(float(mpmath_pi) - my_pi)/float(mpmath_pi)
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
3.141592653589794
pi is approximately: 3.141592653589794 (as a float number, precision is limited). Accuracy % with math.pi: 2.827e-14 Accuracy % with mpmath_pi: 2.827e-14
my_pi = gauss_legendre_1(40)
my_pi
print("pi is approximately: {:.15f} (as a float number, precision is limited).".format(my_pi))
accuracy = 100*abs(math.pi - my_pi)/math.pi
print("Accuracy % with math.pi: {:.4g}".format(accuracy))
accuracy = 100*abs(float(mpmath_pi) - my_pi)/float(mpmath_pi)
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
3.141592653589794
pi is approximately: 3.141592653589794 (as a float number, precision is limited). Accuracy % with math.pi: 2.827e-14 Accuracy % with mpmath_pi: 2.827e-14
This first implementation of the Gauss-Legendre algorithm is limited to a precision of 13 or 14 digits. But it converges quickly ! (4 steps here).
from decimal import Decimal, getcontext
def gauss_legendre_2(max_step):
"""Decimal number implementation of the Gauss-Legendre algorithm, for max_step steps."""
# trick to improve precision
getcontext().prec = 3 + 2**(max_step + 2)
cst_2 = Decimal(2.0)
cst_4 = Decimal(4.0)
a = Decimal(1.0)
b = Decimal(0.5).sqrt()
t = Decimal(0.25)
p = Decimal(1.0)
for i in range(max_step):
new_a = (a+b)/cst_2
new_b = (a*b).sqrt()
new_t = Decimal(t - p*(a - new_a)**2)
new_p = cst_2*p
a, b, t, p = new_a, new_b, new_t, new_p
my_pi = Decimal(((a+b)**2)/(cst_4*t))
return my_pi
my_pi = gauss_legendre_2(5)
print("pi is approximately: {}.".format(my_pi.to_eng_string()[:2**(5+1)]))
accuracy = 100*abs(Decimal(math.pi) - my_pi)/Decimal(math.pi)
print("Accuracy % with math.pi: {:.4g}".format(accuracy))
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
pi is approximately: 3.14159265358979323846264338327950288419716939937510582097494459. Accuracy % with math.pi: 3.898e-15 Accuracy % with mpmath_pi: 7.659e-83
The second implementation of the Gauss-Legendre algorithm is now working better (when we adapt the precision). And it converges quickly ! (8 steps give a precision upto the 697th digits).
How much did we lost in term of performance by using decimal numbers? About a factor $1000$, but that's normal as the second solution stores a lot of digits. They don't even compute the same response.
%timeit gauss_legendre_1(8)
%timeit gauss_legendre_2(8)
The slowest run took 6.47 times longer than the fastest. This could mean that an intermediate result is being cached. 100000 loops, best of 3: 4.06 µs per loop 100 loops, best of 3: 3.9 ms per loop
Methods based on a convergent series¶
For the following formulae, you can try to write a program that computes the partial sum of the series, up to a certain number of term ($N \geq 1$). Basically, the bigger the $N$, the more digits you get (but the longer the program will run).
Some methods might be really efficient, only needing a few number of steps (a small $N$) for computing many digits.
Try to implement and compare at least two of these methods.
You can use the function sum (or math.fsum) to compute the sum, or a simple while/for loop.
All these partial sums are written up to an integer $N \geq 1$.
A Leibniz formula (easy):¶
It has a number of digits sub-linear in the number $N$ of terms in the sum: we need $10$ times more terms to win just one digit. $$\pi \simeq 4\sum_{n=0}^{N} \cfrac {(-1)^n}{2n+1}. $$
def leibniz(max_step):
""" Computing an approximation of pi with Leibniz series."""
my_pi = Decimal(0)
for k in range(max_step):
my_pi += Decimal((-1)**k) / Decimal(2*k+1)
return Decimal(4) * my_pi
getcontext().prec = 20 # trick to improve precision
my_pi = leibniz(1000)
my_pi
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
Decimal('3.1405926538397929257')
Accuracy % with mpmath_pi: 0.03183
getcontext().prec = 20 # trick to improve precision
my_pi = leibniz(10000)
my_pi
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
Decimal('3.1414926535900432389')
Accuracy % with mpmath_pi: 0.003183
This first formula is very inefficient!
Bailey-Borwein-Plouffe series (medium):¶
It also has a number of digits linear in the number $N$ of terms in the sum. $$\pi \simeq \sum_{n = 1}^{N} \left( \frac{1}{16^{n}} \left( \frac{4}{8n+1} - \frac{2}{8n+4} - \frac{1}{8n+5} - \frac{1}{8n+6} \right) \right). $$
def bbp(max_step):
""" Computing an approximation of pi with Bailey-Borwein-Plouffe series."""
my_pi = Decimal(0)
for k in range(max_step):
my_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)))
return my_pi
getcontext().prec = 20 # trick to improve precision
my_pi = bbp(10)
my_pi
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
Decimal('3.1415926535897911465')
Accuracy % with mpmath_pi: 6.659e-14
That's pretty impressive, in only $10$ steps! But that algorithm is highly dependent on the precision we ask, and on the number of terms in the sum.
getcontext().prec = 200 # trick to improve precision
my_pi = bbp(200)
my_pi
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303817')
Accuracy % with mpmath_pi: 8.417e-198
getcontext().prec = 500 # trick to improve precision
my_pi = bbp(500)
my_pi
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119494')
Accuracy % with mpmath_pi: 8.600e-498
It is, of course, slower than the optimized algorithm from mpmath.
And it does not scale well with the precision:
getcontext().prec = 10 + 1000 # trick to improve precision
%timeit bbp(1000)
getcontext().prec = 10 + 2000 # trick to improve precision
%timeit bbp(2000)
10 loops, best of 3: 95 ms per loop 1 loop, best of 3: 661 ms per loop
Bellard's formula (hard):¶
It is a more efficient formula. $$\pi \simeq \frac1{2^6} \sum_{n=0}^N \frac{(-1)^n}{2^{10n}} \, \left(-\frac{2^5}{4n+1} - \frac1{4n+3} + \frac{2^8}{10n+1} - \frac{2^6}{10n+3} - \frac{2^2}{10n+5} - \frac{2^2}{10n+7} + \frac1{10n+9} \right). $$
def bellard(max_step):
""" Computing an approximation of pi with Bellard series."""
my_pi = Decimal(0)
for k in range(max_step):
my_pi += (Decimal(-1)**k/(1024**k))*( Decimal(256)/(10*k+1) + Decimal(1)/(10*k+9) - Decimal(64)/(10*k+3) - Decimal(32)/(4*k+1) - Decimal(4)/(10*k+5) - Decimal(4)/(10*k+7) -Decimal(1)/(4*k+3))
return my_pi * Decimal(1.0/(2**6))
getcontext().prec = 40 # trick to improve precision
my_pi = bellard(10)
my_pi
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
Decimal('3.141592653589793238462643383279490036578')
Accuracy % with mpmath_pi: 4.090e-31
That's pretty impressive, in only $10$ steps! But that algorithm is again highly dependent on the precision we ask, and on the number of terms in the sum.
getcontext().prec = 800 # trick to improve precision
my_pi = bellard(200)
my_pi
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005611526979981607242018469874601593584399987639992795109821175268982496078980620087198769557504210588066903133631260918929384314932145237223424648852944318026933908108945450734428224399103245946394')
Accuracy % with mpmath_pi: 2.220e-604
It is, of course, slower than the optimized algorithm from mpmath.
And it does not scale well with the precision:
getcontext().prec = 10 + 1000 # trick to improve precision
%timeit bellard(1000)
getcontext().prec = 10 + 2000 # trick to improve precision
%timeit bellard(2000)
1 loop, best of 3: 378 ms per loop 1 loop, best of 3: 2.91 s per loop
It is also slower than BBP formula.
One Ramanujan's formula (hard):¶
It is efficient but harder to compute easily with high precision, due to the factorial.
But hopefully, the function math.factorial works with Python integers, of arbitrary size, so this won't be an issue.
Remark: This method and the next one compute approximation of $\frac{1}{\pi}$, not $\pi$.
This method has a quadratic precision (the number of correct digits is of the order $\mathcal{O}(N^2)$.
from math import factorial
def ramanujan(max_step):
""" Computing an approximation of pi with a Ramanujan's formula."""
my_pi = Decimal(0)
d_1103 = Decimal(1103)
d_26390 = Decimal(26390)
d_396 = Decimal(396)
for k in range(max_step):
my_pi += ((Decimal(factorial(4 * k))) * (d_1103 + d_26390 * Decimal(k))) / ( (Decimal(factorial(k)))**4 * (d_396**(4*k)))
my_pi = my_pi * 2 * Decimal(2).sqrt() / Decimal(9801)
my_pi = my_pi**(-1)
return my_pi
getcontext().prec = 40 # trick to improve precision
my_pi = ramanujan(4)
my_pi
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
Decimal('3.141592653589793238462643383279555273157')
Accuracy % with mpmath_pi: 1.668e-30
getcontext().prec = 400 # trick to improve precision
my_pi = ramanujan(40)
my_pi
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
Decimal('3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588249728368707403112841783926380916831222457885852997711270901024127742895829505936')
Accuracy % with mpmath_pi: 2.382e-318
getcontext().prec = 2000 # trick to improve precision
my_pi = ramanujan(200)
my_pi
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989380952572010654858632788659361533818279682303019520353018529689957736225994138912497217752834791315155748572424541506959508295331168617278558890750983817546374649393192550604009277016711390098488240128583616035637076601047101819429555961989467678374494482553797747268471040475346462080466842590694912933136770289891521047521620569660240580381501935112533824300355876402474964732639141992726042699227967823547816360093417216412199245863150302861829745557067498385054945885869269956909272107975093029553211653449872027559602364806654991198818347977535663698074265425278625518184175746728909777727940092593679577412880114072520662452344498309284389026310583057432459663517241665671357564259409181627161036619398380266804323123256475817186916514676524019320976885611753635255011894984447618075599914209997043044425803288826839697424724950586712681051712136647565220069328838660790206792014356556715545181106432983263815336472301873472006826126207435983354628153181952126467499415321188698738323941189429')
Accuracy % with mpmath_pi: 6.658e-1596
$1595$ correct digits with $200$ terms, that's quite good!!
getcontext().prec = 10 + 2000 # trick to improve precision
%timeit ramanujan(200)
getcontext().prec = 10 + 5000 # trick to improve precision
%timeit ramanujan(400)
10 loops, best of 3: 50.1 ms per loop 1 loop, best of 3: 478 ms per loop
Let's try to answer my initial question, using this naive implementation.
%%time
getcontext().prec = 140350 # trick to improve precision
i = 140317
my_pi = ramanujan(10000)
print(str(my_pi)[2:][i:i+10])
mp.dps=140330
print(str(mp.pi)[2:][i:i+10])
--------------------------------------------------------------------------- KeyboardInterrupt Traceback (most recent call last) <ipython-input-128-f3c69c467c77> in <module>() ----> 1 get_ipython().run_cell_magic('time', '', 'getcontext().prec = 140350 # trick to improve precision\ni = 140317\nmy_pi = ramanujan(10000)\nprint(str(my_pi)[2:][i:i+10])\n\nmp.dps=140330\nprint(str(mp.pi)[2:][i:i+10])') /usr/local/lib/python3.5/dist-packages/IPython/core/interactiveshell.py in run_cell_magic(self, magic_name, line, cell) 2113 magic_arg_s = self.var_expand(line, stack_depth) 2114 with self.builtin_trap: -> 2115 result = fn(magic_arg_s, cell) 2116 return result 2117 <decorator-gen-59> in time(self, line, cell, local_ns) /usr/local/lib/python3.5/dist-packages/IPython/core/magic.py in <lambda>(f, *a, **k) 186 # but it's overkill for just that one bit of state. 187 def magic_deco(arg): --> 188 call = lambda f, *a, **k: f(*a, **k) 189 190 if callable(arg): /usr/local/lib/python3.5/dist-packages/IPython/core/magics/execution.py in time(self, line, cell, local_ns) 1183 else: 1184 st = clock2() -> 1185 exec(code, glob, local_ns) 1186 end = clock2() 1187 out = None <timed exec> in <module>() <ipython-input-123-e46be6bb0995> in ramanujan(max_step) 8 d_396 = Decimal(396) 9 for k in range(max_step): ---> 10 my_pi += ((Decimal(factorial(4 * k))) * (d_1103 + d_26390 * Decimal(k))) / ( (Decimal(factorial(k)))**4 * (d_396**(4*k))) 11 my_pi = my_pi * 2 * Decimal(2).sqrt() / Decimal(9801) 12 my_pi = my_pi**(-1) KeyboardInterrupt:
... It was too slow!
Chudnovsky brothers' formula (hard):¶
$$\frac{1}{\pi} \simeq 12 \sum_{n=0}^N \frac {(-1)^{n}(6n)!(545140134n+13591409)}{(3n)!(n!)^{3}640320^{{3n+3/2}}}. $$In 2015, the best method is still the Chudnovsky brothers's formula.
Be careful when you use these formulae, check carefully the constants you wrote (545140134 will work well, 545140135 might not work as well!).
from math import factorial
def chudnovsky(max_step):
""" Computing an approximation of pi with Bellard series."""
my_pi = Decimal(0)
for k in range(max_step):
my_pi += (Decimal(-1)**k)*(Decimal(factorial(6*k))/((factorial(k)**3)*(factorial(3*k)))* (13591409+545140134*k)/(640320**(3*k)))
my_pi = my_pi * Decimal(10005).sqrt()/4270934400
my_pi = my_pi**(-1)
return my_pi
It is very efficient, as Ramanujan's formula.
getcontext().prec = 3000 # trick to improve precision
my_pi = chudnovsky(200)
my_pi
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
Decimal('3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778053217122680661300192787661119590921642019893809525720106548586327886593615338182796823030195203530185296899577362259941389124972177528347913151557485724245415069595082953311686172785588907509838175463746493931925506040092770167113900984882401285836160356370766010471018194295559619894676783744944825537977472684710404753464620804668425906949129331367702898915210475216205696602405803815019351125338243003558764024749647326391419927260426992279678235478163600934172164121992458631503028618297455570674983850549458858692699569092721079750930295532116534498720275596023648066549911988183479775356636980742654252786255181841757467289097777279380008164706001614524919217321721477235014144197356854816136115735255213347574184946843852332390739414333454776241686251898356948556209921922218427255025425688767179049460165346680498862723279178608578438382796797668145410095388378636095068006422512520511739298489608412848862694560424196528502221066118630674427862203919494504712371378696095636437191728746776465757396241389086583264599581339047802759009946576407895126946839835259570982582262052248940772671947826848260147699090264013639443745530506820349625245174939965143142980919065925093722169646151570985838741059788595977297549893016175392846813826868386894277415599185592524595395943104997252468084598727364469584865383673622262609912460805124388439045124413654976278079771569143599770012961608944169486855584840635342207222582848864815845602850601684273945226746767889525213852254995466672782398645659611635488623057745649803559363456817432411251507606947945109659609402522887971089314566913686722874894056010150330861792868092087476091782493858900971490967598526136554978189312978482168299894872265880485756401427047755513237964145152374623436454285844479526586782105114135473573952311342716610213596953623144295248493718711014576540359027993440374200731057853906219838744780847852710405603214833731504646600861733152271806431614075966417958570002328166126409832551776784733662072536424251339850953855478119486626262729703543651855423452347092')
Accuracy % with mpmath_pi: 1.191e-2835
It gets $2834$ correct numbers in $200$ steps! It is more efficient that Ramanujan's formula, but slower to compute.
getcontext().prec = 6000 # trick to improve precision
my_pi = chudnovsky(400)
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
Accuracy % with mpmath_pi: 3.947e-5672
getcontext().prec = 3000 # trick to improve precision
%timeit chudnovsky(200)
getcontext().prec = 6000 # trick to improve precision
%timeit chudnovsky(400)
1 loop, best of 3: 216 ms per loop 1 loop, best of 3: 1.93 s per loop
About $2$ seconds to find correctly the first $5671$ digits? That's slow! But hey, it's Python (dynamic typing etc).
Other methods¶
Trigonometric methods (hard)¶
Some methods are based on the $\mathrm{arccot}$ or $\arctan$ functions, and use the appropriate Taylor series to approximate these functions. The best example is Machin's formula: $$\pi = 16 \;\mathrm{arccot}(5) - 4 \;\mathrm{arccot}(239).$$
And we use the Taylor series to approximate $\mathrm{arccot}(x)$: $$\mathrm{arccot}(x) = \frac{1}{x} - \frac{1}{3x^3} + \frac{1}{5x^5} - \frac{1}{7x^7} + \dots = \sum_{n=0}^{+\infty} \frac{(-1)^n}{(2n+1) x^{2n+1}} .$$
This method is also explained here with some details.
In order to obtain $n$ digits, we will use fixed-point arithmetic to compute $\pi \times 10^n$ as a Python long integer.
High-precision arccot computation¶
To calculate $\mathrm{arccot}$ of an argument $x$, we start by dividing the number $1$ (represented by $10^n$, which we provide as the argument unity) by $x$ to obtain the first term.
We then repeatedly divide by $x^2$ and a counter value that runs over $3$, $5$, $7$ etc, to obtain each next term.
The summation is stopped at the first zero term, which in this fixed-point representation corresponds to a real value less than $10^{-n}$:
def arccot(x, unity):
xpower = unity / x
sum = xpower
n = 3
sign = -1
while True:
xpower = xpower / (x*x)
term = xpower / n
if term == 0:
break # we are done
sum += sign * term
sign = -sign
n += 2
return sum
Adapting it to use Decimal numbers is easy:
def arccot(x, unity):
"""Compute arccot(x) with a certain level of precision."""
x = Decimal(x)
unity = Decimal(unity)
mysum = xpower = unity / x
n = 3
sign = -1
while True:
xpower = xpower / (x*x)
term = xpower / n
if not term:
break
mysum += sign * term
sign = -sign # we alternate the sign
n += 2
return mysum
Applying Machin's formula¶
Finally, the main function uses Machin's formula to compute $\pi$ using the necessary level of precision, by using this high precision $\mathrm{arccot}$ function: $$\pi = 16 \;\mathrm{arccot}(5) - 4 \;\mathrm{arccot}(239).$$
def machin(digits):
unity = 10**(digits + 10)
pi = 4 * (4*arccot(5, unity) - arccot(239, unity))
return pi / unity
To avoid rounding errors in the result, we use 10 guard digits internally during the calculation. We may now reproduce the historical result obtained by Machin in 1706.
def machin(digits):
"""Compute pi with Machin's formula, with precision at least digits."""
unity = 10**(digits + 10)
my_pi = Decimal(4) * (Decimal(4)*arccot(5, unity) - arccot(239, unity))
return my_pi / Decimal(unity)
getcontext().prec = 10000 # trick to improve precision
my_pi = machin(100)
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
Accuracy % with mpmath_pi: 1.501e-9996
So we got the first $9995$ digits correctly... in $45$ seconds. That will still be too slow to have the digits at position $130317$.
getcontext().prec = 5000 # trick to improve precision
%timeit machin(50)
getcontext().prec = 10000 # trick to improve precision
%timeit machin(100)
1 loop, best of 3: 22.7 s per loop 1 loop, best of 3: 44.5 s per loop
The program can be used to compute tens of thousands of digits in just a few seconds on a modern computer. Typical implementation will be in highly efficient compiled language; like C or maybe Julia.
Many Machin-like formulas also exist, like: $$\pi = 4\arctan\left(\frac{1}{2}\right) + 4 \arctan\left(\frac{1}{3}\right).$$
Trying to solve my question!¶
The important parameter to tune is not the precision given to machin() but the Decimal.prec precision.
%%time
i = 14031
getcontext().prec = i + 20 # trick to improve precision
mp.dps = i + 20
print(str(mp.pi)[2:][i:i+10])
my_pi = machin(11)
print(str(my_pi)[2:][i:i+10])
7126218283 7126218283 CPU times: user 1min 3s, sys: 12 ms, total: 1min 3s Wall time: 1min 3s
%%time
i = 140317
getcontext().prec = i + 20 # trick to improve precision
mp.dps = i + 20
print(str(mp.pi)[2:][i:i+10])
my_pi = machin(50)
print(str(my_pi)[2:][i:i+10])
9341076406 9341076406 CPU times: user 10min 40s, sys: 224 ms, total: 10min 40s Wall time: 10min 41s
It was too slow too! But at least it worked!
My manual implementation of Machin's formula took about $10$ minutes, on an old laptop with $1$ core running Python 3.5, to find the $10$ digits of $\pi$ at index $140317$.
Conclusion¶
$\implies$ To the question "What are the $10$ digits of $\pi$ at index $140317..140326$", the answer is $9341076406$.
(hard) Unbounded Spigot Algorithm¶
This algorithm is quite efficient, but not easy to explain. Go check on-line if you want.
See this page (http://codepad.org/3yDnw0n9) for a 1-line C program that uses a simpler Spigot algorithm for computing the first 15000 digits
A nice method, with a generator yielding the next digit each time, comes from http://stackoverflow.com/a/9005163.
It uses only integer, so no need to do any Decimal trick here.
def next_pi_digit(max_step):
q, r, t, k, m, x = 1, 0, 1, 1, 3, 3
for j in range(max_step):
if 4 * q + r - t < m * t:
yield m
# More details on Python generators can be found here http://stackoverflow.com/a/231855
q, r, t, k, m, x = 10*q, 10*(r-m*t), t, k, (10*(3*q+r))//t - 10*m, x
else:
q, r, t, k, m, x = q*k, (2*q+r)*x, t*x, k+1, (q*(7*k+2)+r*x)//(t*x), x+2
def generator_pi(max_step):
big_str = ''.join(str(d) for d in next_pi_digit(max_step))
return Decimal(big_str[0] + '.' + big_str[1:])
It does not use Decimal numbers.
getcontext().prec = 50 # trick to improve precision
generator_pi(1000)
Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337')
getcontext().prec = 5000 # trick to improve precision
generator_pi(1000)
Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337')
(hard) Borwein's algorithm¶
It has several versions, one with a cubic convergence (each new step multiplies by $3$ the number of digits), one with a nonic convergence (of order $9$) etc. They are not so easy to explain either. Please check on-line, here en.WikiPedia.org/wiki/Borwein's_algorithm.
The cubic method is similar to the Gauss-Legendre algorithm:
- Start with $a_0 = 1/3$, and $s_0 = \frac{\sqrt{3}-1}{2}$,
- And then iterate, as much as you want, by defining $r_{k+1} = \frac{3}{1+2(1-s_k^3)^{1/3}}$, and updating $s_{k+1} = \frac{r_{k+1}-1}{2}$ and $a_{k+1} = r_{k+1}^2 a_k - 3^k (r_{k+1}^2 - 1)$.
Then the numbers $a_k$ will converge to $\frac{1}{\pi}$.
Examples and references¶
Links¶
- en.WikiPedia.org/wiki/Pi#Modern_quest_for_more_digits,
- www.JoyOfPi.com/pi.html and www.JoyOfPi.com/pilinks.html,
- www.EveAndersson.com/pi/digits/ has great interactive tools,
- more crazy stuff MathWorld.Wolfram.com/PiDigits.html, or MathWorld.Wolfram.com/Pi.html,
- one idea with Fibonacci numbers,
- and this incredibly long list of digits at PiWorld.calico.jp/estart.html.
- http://bellard.org/pi/ by Francois Bellard (and http://bellard.org/pi/pi_n2/pi_n2.html)
That's it for today! Happy Pi Day!
