Approximating $\pi$ using Euler's identity with primes
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]):
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
Enjoy Reading This Article?
Here are some more articles you might like to read next: