The Basel problem and Euler’s product formula

Euler showed in 1735 that \(\sum 1/n^2 = \pi^2/6\), and in 1737 that the same sum is a product over the primes, \(\zeta(s) = \prod_p (1 - p^{-s})^{-1}\). This script compares how fast the partial sums and the truncated Euler product approach \(\pi^2/6\), and shows the product diverging at \(s = 1\).

import math

import matplotlib.pyplot as plt
import numpy as np
from scipy.special import zeta

from mathematicskit.number_theory import euler_product

Partial sums vs. the Euler product at s = 2

limits = np.unique(np.logspace(1, 6, 30).astype(int))
partial_sums = np.array([np.sum(1.0 / np.arange(1, n + 1) ** 2) for n in limits])
products = np.array([euler_product(2.0, int(n)) for n in limits])
print(f"pi^2/6               = {math.pi**2 / 6:.10f}")
print(f"scipy zeta(2)        = {zeta(2.0):.10f}")
print(f"Euler product p<=1e6 = {products[-1]:.10f}")

fig, ax = plt.subplots()
ax.loglog(limits, math.pi**2 / 6 - partial_sums, "o-", label="sum over n <= N")
ax.loglog(limits, np.abs(math.pi**2 / 6 - products), "s-", label="product over primes p <= N")
ax.set_xlabel("N")
ax.set_ylabel("|error| vs. pi^2 / 6")
ax.set_title("Two roads to zeta(2)")
ax.legend()

pi^2/6               = 1.6449340668
scipy zeta(2)        = 1.6449340668
Euler product p<=1e6 = 1.6449339554

<matplotlib.legend.Legend object at 0x11a2a2510>

At s = 1 the product diverges: infinitely many primes

Mertens’ third theorem gives \(\prod_{p \le N}(1-1/p)^{-1} \sim e^{\gamma}\ln N\); a finite set of primes would give a finite product.

harmonic = np.array([euler_product(1.0, int(n)) for n in limits])
fig, ax = plt.subplots()
ax.semilogx(limits, harmonic, "o-", label="product over p <= N at s = 1")
ax.semilogx(limits, np.exp(np.euler_gamma) * np.log(limits), "--", color="firebrick", label="e^gamma ln N")
ax.set_xlabel("N")
ax.legend()
ax.set_title("Euler's proof that the primes are infinite")
plt.show()