π approximation with Python and mpmath

Computers were originally made for computing, hence the name - computers. In this post we are going to try to compute the π constant with at least medium precision, let's say, to the first 150 digits.

There is a wide choice of π approximation methods: https://en.wikipedia.org/wiki/Approximations_of_%CF%80

People made attempts to approximate the number π since Before Christ, so the approaches range from ancient to really modern. We are going to use a method devised by a French mathematician François Viète in 16th century: https://en.wikipedia.org/wiki/Vi%C3%A8te%27s_formula

We are going to use Python to implement the formula. First thing we can notice is, the numerators of each of the multipliers in the formula nicely depend on each other. The 2nd numerator is the square root of 2 + the previous numerator. The 3rd one is the square root of 2 + the 2nd numerator, and so on. It means we could generate an array of the numerators:

def viete_numerators(n):

    current = 1

    result = []

    for i in range(n):

        if i == 0:

            current = sqrt(2)

        else:

            current = sqrt(2 + current)

        result.append(current)

    return result

Then we can add the division by two to each element of the array:

map(lambda x: x/2, viete_numerators(n))

and then use functools.reduce to multiply the elements of the array:

reduce(lambda x,y: x*y, map(lambda x: x/2, viete_numerators(n)))

So, our calculation of π is:

def pi_approx(n):

    multiplication = reduce(lambda x,y: x*y, map(lambda x: x/2, viete_numerators(n)))

    return 2 / multiplication

Let's try running it with n = 10 (i.e. we generate 10 elements of the Viete's formula):

3.1415914215111997 - nice, but the result is correct only up to 3.14159. What if we do the same for n = 1000?

3.141592653589794 - much better! Only the last digit was rounded - in better approximations, it is 3, not four. How do we get better approximations then? Well, increasing the n will not help here anymore, because we are constrained by the precision of the float type in Python. We need another numerical type, better suited to what we are trying to achieve here.

One of the options is to use the mpmath Python module. The module allows us to arbitrarily set the precision.

The complete code (also available at https://github.com/pgorak/pi-approximation)

from functools import reduce

from mpmath import mp

def viete_numerators(n):

    current = 1

    result = []

    for i in range(n):

        if i == 0:

            current = mp.sqrt(2)

        else:

            current = mp.sqrt(2 + current)

        result.append(current)

    return result

def pi_approx(n):

    multiplication = reduce(lambda x,y: x*y, map(lambda x: x/2, viete_numerators(n)))

    return 2 / multiplication

if __name__ == '__main__':

    mp.dps = 150

    print(pi_approx(1000))

Please note the three uses of mpmath:

• Two uses of mp.sqrt - we are not using the regular math.sqrt function, but the mpmath's one
  
• The setting of the precision in "main"

The result of running it is:

3.14159 26535 89793 23846 26433 83279 50288 41971 69399 37510 58209 74944 59230 78164 06286 20899 86280 34825 34211 70679 82148 08651 32823 06647 09384 46095 50582 23172 53594 0813

and this is a decent approximation - only the last digit "3" is incorrect (rounded). You can now run it a few times, trying to decrease the value of n to see how many elements of the Viete's formula are enough to produce the same result (we don't need as many as 1000).

This example shows that for high precision calculations we need something more sophisticated than the native types in a language.