Estimating Pi in Python

For this project I will code a few of the many methods of estimating Pi, the incredibly useful number that tells us the ratio of a circle's circumference to its diameter.

Nov 05, 2024

Pi is an irrational number starting off 3.14159 and then carrying on for an infinite number of digits with no pattern which anybody has ever discovered. Therefore it cannot be calculated as such, just estimated to (in principle) any number of digits.

So far π has been calculated to 22,459,157,718,361 digits (nearly twenty two and a half trillion) which, by coincidence, is the same number of different ways there are of calculating it. OK, that's a slight exaggeration but there are many ways of getting the job done, ranging from the ancient 22/7 to recently discovered fast but complex methods such as the Chudnovsky Algorithm. These generally calculate all the digits up to a certain point but there are also a few so-called spigot algorithms which calculate a given digit without calculating all the preceding ones first.

There is plenty of information around regarding the number of digits of π which are actually necessary for practical purposes. This caught my eye though from NASA. You can calculate the circumference of the Universe to the accuracy of a hydrogen atom using just 40dp which is the number reached over three centuries ago! And you only need 62 digits to calculate the circumference to the Planck length, the smallest meaningful unit of length which you can think of as a quantum of space.

Trillions of digits aren't therefore useful in their own right but the effort put in to their calculation no doubt benefits mathematics and computer science.

The Code

All the code for this project is in a single file called estimatingpi.py which you can clone or download from the Github repository.

This is the first part of the source code.

estimatingpi.py part 1

import math


PI_STRING = "3.141592653589793238"
RED = "\x1B[31m"
GREEN = "\x1B[32m"
RESET = "\x1B[0m"


def main():

    """
    Here we simply call the functions which 
    estimate pi using various methods
    """

    print("-----------------")
    print("| codedrome.com |")
    print("| Estimating Pi |")
    print("-----------------\n")

    fractions()

    #francois_viete()

    #john_wallis()

    #john_machin()

    #gregory_leibniz()

    #nilakantha()


def print_as_text(pi):

    """
    Takes a value for pi and prints it below a definitive value,
    with matching digits in green and non-matching digits in red
    """

    pi_string = str("%1.18f" % pi)

    print("Definitive: " + PI_STRING)

    print("Estimated:  ", end="")

    for i in range(0, len(pi_string)):

        if pi_string[i] == PI_STRING[i]:

            print(GREEN +  pi_string[i] + RESET, end="")

        else:

            print(RED +  pi_string[i] + RESET, end="")

    print("\n")

Firstly we import math. Not a big surprise! Then I have created a string to hold Pi to 18dp; this is used to compare digit by digit with our calculated values. Next we have a few ANSI terminal codes for printing in colour.

Within main we have calls to all the functions which we'll implement. We can uncomment them one at a time as we write the functions.

Next comes a quick utility function called print_as_text before we get into actually estimating pi, which takes a value for pi, uses it to create a string to 18dp, and prints it out below the definitive value for comparison, with matching digits in green and non-matching ones in red.

Estimating Pi Using Fractions

The first of the six methods we will use is the simplest and just consists of dividing one number by another, or to look at it another way converting fractions to decimals

estimatingpi.py part 2

def fractions():

    """
    Estimates pi using a selection of increasingly accurate fractions
    """

    pi = 22 / 7
    print("22/7\n====")
    print_as_text(pi)

    pi = 333 / 106
    print("333/106\n=======")
    print_as_text(pi)

    pi = 355 / 113
    print("355/113\n=======")
    print_as_text(pi)

    pi = 52163 / 16604
    print("52163/16604\n===========")
    print_as_text(pi)

    pi = 103993 / 33102
    print("103993/33102\n============")
    print_as_text(pi)

    pi = 245850922 / 78256779
    print("245850922/78256779\n==================")
    print_as_text(pi)

The variable pi is set to each of six different fractions and then printed using the print_as_text function. Uncomment fractions in main and run the program using this command...

python3 estimatingpi.py...

to get this.

As you can see the longer fractions give us ever increasing accuracy. The second fraction looks odd as it has a seemingly correct digit after a number of wrong ones but that it just a coincidence. It is quite impressive that pi can be calculated to such accuracy with a simple long division, although I suspect the later fractions were discovered after pi had already been calculated to that accuracy using more complex methods.

Francois Viete

The following formula was discovered by French mathematician Francois Viete in 1593.

These are just the first three terms of an infinite series and already grow to somewhat monstrous proportions. However, if you look closely you can see that each numerator (the expression on the top of the fraction) is just the square root of 2 plus the previous numerator. This even applies to the first term if you assume the previous numerator is 0, so let's use that insight to implement Monsieur Viete's infinite product.

estimatingpi.py part 3

def francois_viete():

    """
    Infinite series discovered by French mathematician 
    Francois Viete in 1593
    """

    print("Francois Viete\n==============")

    iterations = 28
    numerator = 0.0
    pi = 1.0

    for i in range(1, iterations + 1):
        numerator = math.sqrt(2.0 + numerator)
        pi*= (numerator / 2.0)

    pi = (1.0 / pi) * 2.0

    print_as_text(pi)

The iterations variable specifies how many terms we will calculate, and the numerator is initialized to 0 as mentioned above. The pi variable is initialised to 1 so we can multiply the value of the first term by it without needing any special case for the first term.

Then we enter a for loop, setting numerator to the square root of 2 + itself as described above, after which pi is multiplied by the evaluated term. Again we do not need a special case for the first term as pi is initialised to 1.

You probably noticed that the formula does not give us π but 2/π so after the loop we need to calculate pi itself and than call print_as_text. Uncomment francois_viete in main and then run again.

Just 28 terms (I guess calculable by hand four hundred years ago with care and patience) give us 14dp of accuracy. Not bad!

John Wallis

English mathematician John Wallis came up with this in 1655, another infinite product (ie. each term is multiplied by the previous), this one giving us π/2 rather than 2/π.

It's easy to see what is going on here, the denominator and numerator are alternately incremented by 2. In the code, if the term is odd we'll increment the denominator, and if it's even we'll increment the numerator.

estimatingpi.py part 4

def john_wallis():

    """
    Infinite product created by English 
    mathematician John Wallis in 1655
    """

    print("John Wallis\n===========")

    iterations = 1000000
    numerator = 2.0
    denominator = 1.0
    pi = 1.0

    for i in range(1, iterations + 1):
        pi*= (numerator / denominator)
        if (i%2) == 1:
            denominator+= 2.0
        else:
            numerator+= 2.0

    pi*= 2.0

    print_as_text(pi)

The iterations variable is set to 1,000,000 (yes, really) and numerator and denominator are set to 2 and 1 respectively, as per the first term. The pi variable is again set to 1 so we can multiply the first term by it without a special case.

The for loop is simple, it just multiplies pi by the next term, incrementing the numerator or denominator depending on whether the current term is odd or even.

So far we only have half of pi so after the loop it's multiplied by 2 and printed. Uncomment john_wallis in main and run.

A million terms, waaaaaay beyond what Mr Wallis could have calculated, and we only get 5dp of accuracy. I described Francois Viete's effort as "not bad" but I have to describe John Wallis's as "not good".

John Machin

This is English astronomer John Machin's contribution from 1706, which he used to calculate 100 digits. Others used his formula to calculate many more digits by hand and it was the most used until the 20th century.

It looks simple - no infinite series, just one calculation. However, the complexity is hidden in the use of arctan although we have Python's math.atan whereas Professor Machin did not.

estimatingpi.py part 5

def john_machin():

    """
    Formula discovered by English astronomer John Machin in 1706
    """

    print("John Machin\n===========")

    pi = (4.0 * math.atan(1.0 / 5.0) - math.atan(1.0 / 239.0)) * 4.0

    print_as_text(pi)

This function is very straightforward, just a translation of the formula into Python. The formula only gives us a measly quarter of π so we need to multiply it by 4 at the end.

Uncomment john_machin in main and run the program.

This looks impressive but of course the accuracy is down to the arctangent we have available. I am not sure how this would have been calculated in Machin's time or how much accuracy he would have achieved.

Gregory-Leibniz

Next we have the Gregory-Leibniz series. This is very simple, it gives us π rather than a fraction or multiple of it and the numerator is constant. All we need to do is add 2 to the denominator each time and flip between subtraction and addition.

estimatingpi.py part 6

def gregory_leibniz():

    """
    Co-discovered by James Gregory and 
    Gottfried Wilhelm Leibniz
    """

    print("Gregory-Leibniz\n===============")

    iterations = 400000
    denominator = 1.0
    multiplier = 1.0
    pi = (4.0 / denominator)

    for i in range(2, iterations + 1):
        denominator += 2.0
        multiplier *= -1.0
        pi += ( (4.0 / denominator) * multiplier )

    print_as_text(pi)

The denominator variable is first set to 1, and we also have a multiplier which will alternate between 1 and -1 to allow for the alternating addition and subtraction. Pi is initialised to the first term, therefore the loop starts at the second term.

Within the loop denominator is incremented by 2 and multiplier is "flipped" before we evaluate and add the next term to pi.

Uncomment the function in main and run again.

I have to admit the result is a bit baffling. The second batch of accurate digits cannot be chance as with 333/106, but why the wildly inaccurate two digits before them?

Nilakantha

Finally another infinite sum, this one from the 15th century Indian mathematician Nilakantha Somayaji, which has some obvious similarities with the Gregory-Leibniz series.

So we start off with 3, have a constant numerator of 4, and calculate the denominator by multiplying three numbers. We also alternate between addition and subtraction.

estimatingpi.py part 7

def nilakantha():

    """
    Named after the 15th century Indian 
    mathematician Nilakantha Somayaji
    """

    print("Nilakantha\n=========")

    iterations = 1000000
    multiplier = 1.0
    start_denominator = 2.0
    pi = 3.0

    for i in range(1, iterations + 1):
        pi += ( (4.0 / (start_denominator * (start_denominator + 1.0) * (start_denominator + 2.0)) ) * multiplier)
        start_denominator += 2.0
        multiplier *= -1.0

    print_as_text(pi)


main()

Again we have a multiplier variable to alternate between addition and subtraction. The start_denominator variable represents the first of the three numbers which are multiplied to get the denominator, and will be incremented by 2 on each iteration. We then start off pi at 3.

Within the loop the numerator 4 is divided by the denominator which is calculated by multiplying start_denominator by its two consective numbers, this then being multiplied by 1 or -1 alternately. Lastly within the loop 2 is added to start_denominator and multiplier is flipped.

Uncomment nilakantha in main and run.

This is pretty good, but even 50 iterations gets us 5dp of accuracy. Remember that each additional digit only gives us at most a 10% increase in accuracy.