Approximating $\pi$ using Euler's identity with primes

dodo · May 21, 2020

Open In Colab

In this notebook, I implement a method to approximate $\pi$ using the identity \(\frac{6}{\pi^2} = \prod_{p} \left(1-\frac{1}{p^2}\right),\) where $p$ are all the prime numbers, starting from 2. If we call the right-hand side of the identity $P$, then we have that \(\pi = \sqrt{\frac{6}{P}}.\)

More information can be found in this YouTube video.

from math import sqrt, pi
import numpy as np

terms = None

The first step of the experiment involves generating a list of all prime numbers up to a number $n$. I got the fastest version of the sieve of Eratosthenes from a bunch of earlier experiments I did.

def eratosthenes(n):
    primes = [2]
    for i in range(3, n, 2):
        if all([i % p != 0 for p in primes]):

    return primes

I created two versions of the function. The first one uses a loop over all primes and accumulates the product in a variable. The second one caches calculated terms of the form $\left(1-\frac{1}{p^2}\right)$ into a numpy array and returns the product of the relevant terms, depending on the function call.

def euler_prime(x):
    a = 1
    for p in eratosthenes(x):
        a *= 1 - (1 / (p ** 2))

    return a

def euler_prime2(x):
    global terms

    num_terms = terms.shape[0] if terms is not None else 0
    primes = eratosthenes(x)
    if num_terms >= len(primes):
        return terms[:len(primes)].cumprod()[-1]

    primes = np.array(primes[num_terms:])
    new_terms = 1 - (1 / (primes ** 2))

    if terms is None:
        terms = new_terms
        terms = np.concatenate((terms, new_terms))

    return terms.cumprod()[-1]
def pi_approx(x):
    return sqrt(6 / euler_prime2(x))

Here are the differences between $\pi$ and out approximations for different numbers of primes.

for i in [3, 10, 183, 2781, 16789]:
    res = abs(pi - pi_approx(i))

Twitter, Facebook