Week 7: Pi

Contents

1 Introduction
2 Computing Pi
3 Efficient Methods
4 Further Reading

This week we look at the irrational number π, and some methods for calculating it.

Introduction

What is Pi?

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²

For a circle, the ratio of the square of its perimeter, P², to its area A, is 4π. For all other shapes, P² / A is greater than 4π. For example, for a square, P² / 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 = 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 and complex numbers

Pi appears in one of the most intriguing equations in mathematics:

Here e = 2.718281... is another irrational number (Euler's number, the base of the natural logarithm), and i is the unit imaginary number (which we encountered in last week's tutorial).

Pi in number theory

π 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 / π² = 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². So the chance of two numbers being coprime is the product of (1 - 1/p²) taken over all primes: i.e., (1 - 1/2²)(1 - 1/3²)(1 - 1/5²) .... Now to simplify this we can use a result from the sum of geometric series:

1 / (1 - 1/p²) = 1 + 1/p² + 1/p⁴ + 1/p⁶ + ....

(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ⁱ q^j r^k ...)²

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² + 1/3² + 1/4² + 1/5² + ... )

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² + 1/3² + 1/4² + 1/5² + ... = π² / 6

A proof of this final step can be found here.

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² + y² < 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)

Exercises

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).

Pi from polygons

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.

Madhava's Formula for π

The introduction of infinite series led to a revolution in calculations of π. The following formula was first discovered by Madhava of Sangamagrama in the 14th century:

π/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:

$\frac{\pi}{4} = 4 \, \arctan \frac{1}{5} - \arctan \frac{1}{239}$

with

$\arctan z = z - \frac {z^3} {3} +\frac {z^5} {5} -\frac {z^7} {7} +\cdots$

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¹⁰⁰ => 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)

in order to find their arithmetic-geometric mean.

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.