[Python] Estimating Pi

 

3.14% of sailors are pi-rates.

 

Pi Estimation

Foto von Chris Linnett auf Unsplash

 

Pi, denoted by the Greek letter π, is an important mathematical constant that represents the ratio of a circle’s circumference to its diameter and is an irrational number, meaning it cannot be expressed as a finite decimal or a fraction. It is usually rounded down to two decimal points as 3.14, however its infinite, non-repeating decimal representation, has fascinated mathematicians and there have been many proposed algorithms for estimating Pi.

Its origins can be traced back to ancient civilizations, where approximations of its value were discovered independently in different cultures. Its symbol, however, was first introduced by the Welsh mathematician William Jones in the 18th century.

Here we will analyze four methods for estimating Pi:

1. Monte Carlo approach

The Monte Carlo approach is an interesting approach for estimating the value of Pi. By utilizing random sampling, this technique provides a clever way to approximate pi’s digits. Although the Monte Carlo approach provides us with a simple way for estimating pi, it still converges rather slowly and is not an optimal solution.

Here you find the code for the Monte Carlo approach:

# Econowmics.com
# Youtube.com/@codedabacus

# libraries
from math import sqrt,pi 
import random 

# The Monte Carlo Approach
def pi_monte_carlo(n):
    # total points in the circle
    c = 0
    
    # generate points
    for _ in range(n):
        x = random.random() #[0,1]
        y = random.random()
        
        # compute the point's distance from the circle
        d = sqrt(x ** 2 + y ** 2)
        
        if d <= 1:
            c += 1
            
    return 4 * (c / n)

print(pi_monte_carlo(1000000))
print(pi)

 

2. Leibniz’s Formula

The Leibniz formula, suggested by the German mathematician Gottfried Wilhelm Leibniz, is another interesting method for approximating the value of pi.  This formula is derived from the Taylor series of the arctangent function. You can read more about its proof here.

Although the convergence is very slow for this approach, it still offers a simple and intuitive approach for estimating Pi. Here you will find the codes for this approach:

# Econowmics.com
# Youtube.com/@codedabacus
from decimal import Decimal, getcontext

# the Leibniz's Formula
def pi_leibniz(K):
    # set the precision
    getcontext().prec = 1000
    
    # starting value
    pi = Decimal('1')
    
    # the summation
    for k in range(1, K):
        # the numerator
        num = Decimal((-1) ** k)
        
        # the denominator
        den = Decimal(2 * k + 1)
        
        pi += Decimal(num / den)
        
    return Decimal(4 * pi)

print(str(pi_leibniz(1000000))[:20])
from math import pi
print(pi)

 

You can find more information about the Monte Carlo approach and the Leibniz’s Formula and their implementation in Python via the video below:

 

 


 

* What you need for the second part:

For the next two approaches I have defined a simple Python function to check for the accuracy of our estimations. Furthermore, I have used a text file, which contains the first actual 1000 digits of Pi, which you can download from here: pi_first_1000.

This is the code for the function for checking the digits of Pi:

# First 1000 digits of Pi
with open(r"pi_first_1000.txt", "r") as f:
    # Read and store the contents of the file
    pi_digits = f.read()

# Function to compare the digits of Pi
def pi_checker(estimate):
    estimate = str(estimate)
    
    for x in range(1000):
        if pi_digits[x] != estimate[x]:
            return x - 2      # because of "3."14...
        
    return 1000

 

3. Bailey-Borwein-Plouffe (BBP) Formula

The Bailey-Borwein-Plouffe (BBP) formula is the third mathematical algorithm which we will analyze for estimating the value of pi. It was named after the three mathematicians who contributed to its discovery: David H. Bailey, Peter Borwein, and Simon Plouffe. The BBP formula provides a novel approach for calculating the hexadecimal digits of pi, which marks it as an extremely useful method for automated calculations. This algorithm converges extremely fast and is very accurate. Find the codes for implementing it here:

# Econowmics.com
# Youtube.com/@codedabacus

from decimal import Decimal, getcontext

# Bailey-Borwein-Plouffe (BBP) Formula
def pi_BBP(K):
    # set the precision
    getcontext().prec = 1000
    
    # starting value
    pi = Decimal(0)
    
    # the summation
    for k in range(K):
        a = Decimal(1) / (16 ** k)
        b = Decimal(4) / (8 * k + 1)
        c = Decimal(2) / (8 * k + 4)
        d = Decimal(1) / (8 * k + 5)
        e = Decimal(1) / (8 * k + 6)
        
        # update the estimate
        pi += Decimal (a * (b - c - d - e))
        
    return pi

res = pi_BBP(1000)
print(str(res)[:17])
from math import pi
print(pi)

print(pi_checker(res))

 

4. Chudnovsky Algorithm

The Chudnovsky algorithm is also a cool method for efficiently estimating the value of pi. Developed by the Chudnovsky brothers, David and Gregory, this algorithm revolutionized the computation of Pi’s decimal places. This algorithm utilizes a series expansion involving factorials, exponentials, and the famous Ramanujan’s formula. With its rapid convergence, the Chudnovsky algorithm allows for the calculation of millions, even billions,  of Pi’s digits. Here you can find the codes for it:

# Econowmics.com
# Youtube.com/@codedabacus

# libraries
from decimal import Decimal, getcontext
from math import factorial

# Chudnovsky Algorithm
def pi_chudnovsky(Q):
    # precision
    getcontext().prec = Q + 50
    
    # starting value
    rhs = 0
    
    # the summation
    for q in range(Q):
        a = (-1) ** q
        b = factorial(6 * q)
        c = (545140134 * q) + 13591409
        d = factorial(3 * q)
        e = factorial(q) ** 3
        f_power = (3 * q) + Decimal(1.5)
        f = 640320 ** f_power
        
        # update our storage
        rhs += Decimal((a * b * c) / (d * e * f))
        
    # return the result
    return 1 / (12 * Decimal(rhs))


res = pi_chudnovsky(1000)
print(pi_checker(res))

 

You can find more information about the BBP Formula and the Chudnovsky Algorithm and their implementation in Python via the video below:

 

Related Images: