Estimate π on π Day

Introduction

Happy π Day 2021! In this article I’ll estimate the digits of π with random numbers and the probability of two integers being co-prime. What is the probability of two random integers being coprime? Euclidean’s Algorithm can be used to estimate π!

Probability of Two Integers Being Coprime

Let denote the sought probability of two random integers being coprime. Pick an integer and two random numbers and . The probability that divides is , and the same holds for . Therefore, the probability that both and are divisible by equals . The probability that and have no other factors, i.e., that equals , by our initial assumption. But is equivalent to . Assuming independence of the events, it follows that the probability that equals .

Now, was just one possibility for the greatest common divisor of two random numbers. Any number could be the . Furthermore, since the events are mutually exclusive (the of two numbers is unique) and the total probability of having a at all is 1 leads to:

implying that

Or denoted in a statistical way

Euclidean Algorithm

In mathematics, the Euclidean Algorithm, or Euclid’s Algorithm, is an efficient method for computing the greatest common divisor (GCD) of two integers (numbers), the largest number that divides them both without a remainder.

Greatest Common Divisor

The greatest common divisor g is the largest natural number that divides both a and b without leaving a remainder.

def gcd(a, b):
    while b != 0:
        t =  b
        b = a % b
        a = t
    return a

Estimate Pi

Let’s step back to the probability of two integers being coprime. We now know that

it can be written as

so we get

Main Code

Talk is cheap. Show me the code.

iterations = 400
num_total = 0
num_coprime = 0
pies = []

for idx in range(iterations):
    for i in range(500):
        a, b = random.randint(1, 1000000), random.randint(1, 1000000)
        if gcd(a, b) == 1:
            num_coprime += 1
        num_total += 1
    pi = math.sqrt(6 * num_total / num_coprime)
    pies.append(pi)
    pi_average = sum(pies)/len(pies)
    print(pi_average)

Visualise

Let’s visualise it dynamicly.

iterations = 400
num_total = 0
num_coprime = 0
pies = []

plt.rcParams["figure.figsize"] = (15, 6)
ax = plt.axes()
plt.xlim(-1, iterations+1) 
plt.ylim(3, 3.3)
plt.axhline(y=3.1415, color='black', linestyle='-')
plt.grid()
for idx in range(iterations):
    for i in range(500):
        a, b = random.randint(1, 1000000), random.randint(1, 1000000)
        if gcd(a, b) == 1:
            num_coprime += 1
        num_total += 1
    pi = math.sqrt(6 * num_total / num_coprime)
    pies.append(pi)
    pi_average = sum(pies)/len(pies)
    ax.plot(
        [idx], [pi_average], 
        linestyle='-', lw=2, marker='o', color='tab:blue', markersize=2, alpha=0.5
    )
    plt.draw()
    custom_lines = [Line2D([0], [0], color='tab:blue', lw=1)]
    plt.legend(custom_lines, [f"pi={pi_average:.10f}"], loc="upper right")
    plt.pause(0.01)
plt.show()

Conclusion

Today I show you how to estimate numbers π from random numbers using Euclidean algorithm. Happy coding! Cheers!

References

1. https://www.cut-the-knot.org/m/Probability/TwoCoprime.shtml
2. https://en.wikipedia.org/wiki/Euclidean_algorithm
3. https://thecodingtrain.com/CodingChallenges/161-pi-from-random.html
4. https://www.youtube.com/watch?v=EvS_a921dBo
5. https://www.youtube.com/watch?v=RZBhSi_PwHU
6. https://www.youtube.com/watch?v=d-o3eB9sfls&list=PLCZeVeoafktVGu9rvM9PHrAdrsUURtLTo&index=54&ab_channel=3Blue1Brown