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.
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:
In [3]:
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))
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).
In [4]:
def mathpi():
from math import pi
return pi
print("First method: using math.pi gives pi ~= {:.17f} (17 digits are displayed here).".format(mathpi()))
If we know the digits, we can directly print them:
In [5]:
from decimal import Decimal
bigpi = Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679')
print("The first 100 digits of pi are {}.".format(bigpi))
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:
In [6]:
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)
In [9]:
print("The simple Monte-Carlo method with {} random points gave pi = {}".format(10000, montecarlo_pi(10000)))
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).
$\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?
mpmathThis 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).
In [10]:
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).
In [11]:
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))
Let save it for further comparison of simpler methods.
In [87]:
mp.dps = 100000 # number of digits
len(str(mp.pi))
mpmath_pi = Decimal(str(mp.pi))
Out[87]:
We can solve the initial challenge easily:
In [86]:
mp.dps = 140330
print(str(mp.pi)[2:][140317:140317+10])
And it will most probably be the quickest method presented here:
In [14]:
%timeit mp.dps=140330;print(str(mp.pi)[2:][140317:140317+10])
Asking for $10$ times more digits take about $100$ more of time (that's a bad news).
In [43]:
%timeit mp.dps=1403230;print(str(mp.pi)[2:][1403217:1403217+10])
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:
Finally, the desired approximation of $\pi$ is: $$\pi \simeq \frac{(a_{n+1} + b_{n+1})^2}{4 t_{n+1}}.$$
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).
In [17]:
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
In [89]:
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))
Out[89]:
In [90]:
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))
Out[90]:
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).
decimal.Decimal to improve precisionThe limitation of this first algorithm came from using simple float numbers.
Let us now use Decimal numbers to keep as many digits after the decimal as we need.
In [24]:
from decimal import Decimal, getcontext
In [25]:
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
In [91]:
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))
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.
In [92]:
%timeit gauss_legendre_1(8)
%timeit gauss_legendre_2(8)
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$.
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}. $$
In [30]:
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
In [98]:
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))
Out[98]:
In [99]:
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))
Out[99]:
This first formula is very inefficient!
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). $$
In [100]:
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
In [101]:
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))
Out[101]:
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.
In [102]:
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))
Out[102]:
In [103]:
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))
Out[103]:
It is, of course, slower than the optimized algorithm from mpmath.
And it does not scale well with the precision:
In [104]:
getcontext().prec = 10 + 1000 # trick to improve precision
%timeit bbp(1000)
getcontext().prec = 10 + 2000 # trick to improve precision
%timeit bbp(2000)
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). $$
In [105]:
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))
In [106]:
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))
Out[106]:
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.
In [107]:
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))
Out[107]:
It is, of course, slower than the optimized algorithm from mpmath.
And it does not scale well with the precision:
In [73]:
getcontext().prec = 10 + 1000 # trick to improve precision
%timeit bellard(1000)
getcontext().prec = 10 + 2000 # trick to improve precision
%timeit bellard(2000)
It is also slower than BBP formula.
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)$.
In [123]:
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
In [124]:
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))
Out[124]:
In [125]:
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))
Out[125]:
In [126]:
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))
Out[126]:
$1595$ correct digits with $200$ terms, that's quite good!!
In [127]:
getcontext().prec = 10 + 2000 # trick to improve precision
%timeit ramanujan(200)
getcontext().prec = 10 + 5000 # trick to improve precision
%timeit ramanujan(400)
Let's try to answer my initial question, using this naive implementation.
In [128]:
%%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])
... It was too slow!
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!).
This formula is used as the logo of the French agrégation of Mathematics website
agreg.org:
In [129]:
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.
In [131]:
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))
Out[131]:
It gets $2834$ correct numbers in $200$ steps! It is more efficient that Ramanujan's formula, but slower to compute.
In [134]:
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))
In [135]:
getcontext().prec = 3000 # trick to improve precision
%timeit chudnovsky(200)
getcontext().prec = 6000 # trick to improve precision
%timeit chudnovsky(400)
About $2$ seconds to find correctly the first $5671$ digits? That's slow! But hey, it's Python (dynamic typing etc).
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.
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:
In [171]:
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
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.
In [172]:
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)
In [173]:
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))
So we got the first $9995$ digits correctly... in $45$ seconds. That will still be too slow to have the digits at position $130317$.
In [174]:
getcontext().prec = 5000 # trick to improve precision
%timeit machin(50)
getcontext().prec = 10000 # trick to improve precision
%timeit machin(100)
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).$$
In [179]:
%%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])
In [180]:
%%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])
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$.
$\implies$ To the question "What are the $10$ digits of $\pi$ at index $140317..140326$", the answer is $9341076406$.
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.
In [149]:
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
In [161]:
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.
In [164]:
getcontext().prec = 50 # trick to improve precision
generator_pi(1000)
Out[164]:
In [165]:
getcontext().prec = 5000 # trick to improve precision
generator_pi(1000)
Out[165]:
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:
Then the numbers $a_k$ will converge to $\frac{1}{\pi}$.
That's it for today! Happy Pi Day!