Approximate Pi with random number generation

approximate-pi.py

# Approximate pi using the probability
# P(gcd(i,j) == 1) ==  6/pi^2
# In other words, the probability that two numbers are co-prime, they have no 
# common denominators, is 6 / (pi^2)
# Or, pi = sqrt(6/P)
#
# Thanks to Matt Parker at Stand Up Maths for this approximation technique, 
# and for awesome content in general: https://www.youtube.com/watch?v=RZBhSi_PwHU
# 
# I didn't feel like making a dedicated repository for this so there is no
# requirements.txt, but it you want to run it please make sure to pip install
# matplotlib and mpmath

import random as r
import math
from mpmath import mpf, mp
from matplotlib import pyplot as plt

def co_prime(n, m):
    return math.gcd(int(n), int(m)) == 1

def approx(p):
    pi = math.sqrt(6 / p)
    return pi

def success(n):
    diff = (mpf(n) - mpf(math.pi)) / mpf(math.pi)
    abs_diff = math.sqrt(diff * diff) # absolute value
    return abs_diff > 0.5

def percent_diff(n, m):
    return (math.fabs(n - m) / m) * 100

def plot(apx):
    successes = [i for (i, n) in enumerate(apx) if math.isclose(n, math.pi, rel_tol=0.0025)]

    plt.plot(apx, label="sqrt(6 / (number of co-primes found))")
    plt.plot([math.pi]*len(apx), label="True pi")
    plt.scatter(successes, [math.pi]*len(successes), label="Percent error < 0.25%", marker="+", c="#51f400")

    plt.title("Approximation of Pi using random number generation")
    plt.ylabel("Approximation")
    plt.xlabel("Random numbers generated")
    plt.legend(loc="best")
    plt.show()

mp.dps = 1000 # 100 bit floating point precision to reduce the number of precision errors, since they get dropped

max_number = 10000

trials = 10000

co_primes = 0

approximations = []

for i in range(trials):
    n = int(r.uniform(1, max_number))
    m = int(r.uniform(1, max_number))

    if co_prime(n, m):
        co_primes = co_primes + 1

    ratio = mpf(co_primes) / mpf(i + 1)
    if ratio == 0 :
        # Precision error. Drop this value
        approximations.append(float('nan'))
    else:
        approximations.append(approx(mpf(ratio)))

print(f"Final approximation of Pi: {approximations[-1]}")
print("Percent error: {}%".format(round(percent_diff(approximations[-1], math.pi), 3)))

plot(approximations)