Pi (π) is a mathematical constant which has played a central role in the history of mathematics:
π = 3.141592653589793238462643383279502884197169 399375105820974944592307816406286208998628034825 342117067982148086513282306647093844609550582231 725359408128481117450284102701938521105559644622 ...
Pi and the circle
π is the ratio of the circumference of a circle (
C) to its diameter (
d):
π = C / d
It's also the ratio of the area of a circle (A) to the square of its radius (r = d/2):
π = A / r^{2}
For a circle, the ratio of the square of its perimeter, P^{2}, to its area A, is 4π. For all other shapes, P^{2} / A is greater than 4π. For example, for a square, P^{2} / A = 16.
Pi is irrational and transcendental
Through the ages, scholars have used simple rational approximations for pi, for example 22/7 or
25/8 or 3927/1250 or 355/113. But π it is
not equal to any fraction: it is an
irrational number. This means that it has an
infinite number of decimal digits, and the digits do not have a repeating pattern.
The square root of 2 is also an irrational number. But unlike π, √2 is a solution to a polynomial equation with integer coefficients:
x^{2} - 2 = 0. For this reason, √2 is called an
algebraic number.
By contrast, Pi is a
transcendental number: it is not the solution to a polynomial equation with integer coefficients.
Pi appears in one of the most intriguing equations in mathematics:
π turned up rather unexpectedly in the the answer to
Mini-Challenge #3 on Farey series. It is intriguing that π makes an appearance here, because Farey series are built from whole numbers, which seem to have little in common with circles and geometry.
Here's a
remarkable fact: the probability of two (large) numbers chosen at random being
relatively prime (or
coprime, that is, having no shared factors) is equal to 6 /
π^{2} = 0.6079... The proof of this statement takes us one step along the path towards one of the most famous unsolved problems in mathematics: the
Riemann hypothesis. Since the steps in the proof can be understood with A-level mathematics, I have included a sketch below:
Proof: Two numbers are not coprime if there exists at least one prime which divides both numbers. The chance of a number being divisible by a given prime
p is 1/
p, and so the chance of both numbers not being divisible by that prime is
1 - 1/p^{2}. So the chance of two numbers being coprime is the product of (1 - 1/p
^{2})
taken over all primes: i.e., (1 - 1/2
^{2})(1 - 1/3
^{2})(1 - 1/5
^{2}) .... Now to simplify this we can use a result from the sum of
geometric series:
1 / (1 - 1/p^{2}) = 1 + 1/p^{2} + 1/p^{4} + 1/p^{6} + ....
(Alternatively, we could use the
binomial theorem). Therefore, each term in the infinite product can be written as (1 divided by) an infinite sum, and, when we multiply out the brackets, we will obtain (in the denominator) a sum of
every possible term of the form:
1 / (p^{i} q^{j} r^{k} ...)^{2}
where p,q,r, etc. are (distinct) primes and i, j, k, etc., are integers greater than or equal to zero. By the
fundamental theorem of arithmetic, every integer greater than 1 has a
unique prime factorization. So every integer will feature in this series. In other words, the chance of two numbers being coprime is:
1 / (1 + 1/2
^{2} + 1/3
^{2} + 1/4
^{2} + 1/5
^{2} + ... )
The problem of finding a
simple solution for the infinite series in the denominator was a famous problem in mathematics called the
Basel problem. It was solved by Leonhard Euler in 1735, who showed that:
1 + 1/2
^{2} + 1/3
^{2} + 1/4
^{2} + 1/5
^{2} + ... =
π^{2} / 6
Computing Pi
A Monte Carlo Method
Let's try calculating π with a simple (but inefficient) method. We will generate some "random" points within the unit square, as shown:
The red points are within the quarter-circle. The blue points are outside the quarter-circle. But how do we tell whether a point with coordinates
x and
y is within the circle? We can use
Pythagoras' theorem to find the distance of the point from the origin. The point is inside the circle
if
x^{2} + y^{2} < 1.
The area of the quarter-circle inside the unit square is π/4. So we expect the proportion of points within the circle (i.e. points shown in red) to be approximately equal to π/4 = 0.785...
This is an example of a
Monte Carlo method.
Such methods rely on repeated random sampling to obtain a numerical
result. It is not a very efficient method. The error in the estimate
decreases only slowly as we increase N, the number of points.
Let's try some Python code:
import random as R
N = 1000
k = 0
count = 0
while k < N:
x, y = R.random(), R.random()
if x**2 + y**2 < 1:
count += 1
k += 1
print "Estimate of Pi using " + str(N) + " points:"
print 4.0*count / float(N)
This code generates N = 1000 pairs of coordinates in the interval 0 < x < 1, 0 < y < 1. It counts how many of these have a length of less than one, and uses this to estimate π.
To generate the coordinates, I'm using the
random module, which implements various (pseudo-)random number generators. Here, the
random()
function returns the next random floating point number in the range [0.0, 1.0).
Of course, with Python you can always write a more condensed program:
import random as R
N = 1000
print 4.0*sum([1 for k in range(N) if R.random()**2 + R.random()**2 < 1])/float(N)
- Show that the error in the estimate decreases with 1/√N, where N is the number of sampling points.
- Write a Monte Carlo program to estimate the proportion of pairs of numbers which are relatively prime, and hence obtain an estimate for π. (First, read the "Pi in number theory" section, above).
Around 250BC, the Greek mathematician Archimedes devised a geometrical approach for calculating Pi, using polygons. π can be estimated by computing the perimeters of circumscribed and inscribed polygons, as below:
Archimedes drew a regular hexagon inside and outside a circle, as above, and then successively
doubled the number of sides until he reached a 96-sided regular
polygon. By calculating the perimeters of the polygons, he proved that
223/71 < π < 22/7.
π/4 = 1 - 1/3 + 1/5 - 1/7 + 1/9 - ...
It is better known as the Leibniz formula for π (alternatively, the Gregory-Leibniz formula).
To use the formula, we must truncate the series. The accuracy of the estimate depends strongly on the rate of convergence of the series. This series converges rather slowly. An approximation for the
error in the estimate is given by the magnitude of the next term that
has been omitted. So if we were to sum the first n terms, the error estimate would be 1/|2n+1|.
We can try Madhava's formula with a one-line Python program:
4.0*sum([1/float(2*k+1)*(-1)**k for k in range(1000)])
Brouncker's Continued Fraction
A rather beautiful formula for Pi was found by Lord William Brouncker in the seventeenth century.
This expression is known as a
generalized continued fraction. Like the Leibniz formula, this formula is slowly converging, and so it is regarded as of little practical use. Yet we may still try it out in Python:
N = 1000
tmp = (2*N+1)**2 / 2.0
k = N - 1
while k >= 0:
tmp = (2*k+1)**2 / (2.0 + tmp)
k -= 1
print 4.0/(1.0 + tmp)
Try increasing N, the number of terms in the continued fraction: how does the error in the estimate change?
Efficient Methods
Now let's look at two efficient methods, which allow us to compute many digits of pi. To compute pi to very high precision, we will need to make use of a Python module: decimal
.
Using the decimal package
The
decimal module provides support for decimal floating point
arithmetic. Crucially, it allows the us to set the
precision of the calculations. That is, we may specify the number of decimal places of accuracy that we need. To see roughly how it works in practice, try the following code:
from decimal import *
setcontext(Context(prec=100))
# Use 100 digits of precision
D = Decimal
D(1.0/3.0)
# using "float"
D(1)/D(3) # using the decimal package
D(2**0.5)
D(2).sqrt()
The former statements, using the standard "float" datatype, should be less accurate than the to the latter statements, using the "decimal" module.
Machin's formula
In 1706
John Machin used the Gregory–Leibniz approach to produce an algorithm that converged much faster than the original:
with
Machin reached 100 digits of
π with this formula. Can we write a Python program to match Machin's achievement? Try this:
import decimal
D = decimal.Decimal
decimal.setcontext(decimal.Context(prec=100))
def arctangent(x, N):
k = 0
tot = D(0)
while k < N:
j = 2*k + 1
tot += (-1)**k * x**j / D(j)
k += 1
return tot
N = 100
print D(4)*( D(4)*arctangent(D(1)/D(5), N) - arctangent(D(1)/D(239), N))
This should give an estimate of pi accurate to 100 decimal places. (When I tested it, all but the last digit seem to be accurate.)
Machin's method is efficient because each term in the series for arctan(1/5) is a factor of at least 25 times smaller than the last. So for 100 digits of accuracy we need only k terms, where 25^{k} = 10^{100} => k = 100 log(10)/log(25) ~ 72.
The Gauss-Legendre Algorithm
An even more efficient method is the
Gauss–Legendre algorithm. It is notable for being rapidly convergent, with only 25 iterations producing
45 million correct digits of π.
It repeatedly replaces two numbers (a & b, say) by their
arithmetic and
geometric means,
(a + b)/2 and √(ab)
The algorithm is quite simple, and is
described here. A Python code for implementing the algorithm, making use of the decimal module, can be
found on Stack Exchange (see the top-rated answer).
I have attached a copy of the code for computing the first 1000 digits of pi at the bottom of this page (pi_gauss_legendre.py
). Please try it for yourself.
Further Reading