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

dodo · May 21, 2020

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}}.$$

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]):
primes.append(i)

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
else:
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))
print(res)

0.3131655288436028
0.04800048589864758
0.0013397623426119054
6.314044652233619e-05
8.757953707139166e-06