Chapter 028 — Prime Numbers

# --- Visualization: Sieve of Eratosthenes and prime density ---
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-whitegrid')

def sieve_eratosthenes(n):
    """Return list of primes up to n using the Sieve of Eratosthenes."""
    is_prime = np.ones(n + 1, dtype=bool)
    is_prime[0] = is_prime[1] = False
    for i in range(2, int(n**0.5) + 1):
        if is_prime[i]:
            is_prime[i*i::i] = False
    return np.where(is_prime)[0]

N = 200
primes = sieve_eratosthenes(N)

fig, axes = plt.subplots(1, 2, figsize=(13, 5))

# Plot 1: Prime grid visualization
ax = axes[0]
grid = np.zeros((14, 15))
for p in primes:
    if p < 14 * 15:
        grid[p // 15, p % 15] = 1
im = ax.imshow(grid, cmap='YlOrRd', aspect='auto', interpolation='nearest')
ax.set_title('Primes up to 210 (dark = prime)')
ax.set_xlabel('Column (mod 15)')
ax.set_ylabel('Row')
for p in primes[:30]:
    r, c = p // 15, p % 15
    ax.text(c, r, str(p), ha='center', va='center', fontsize=7, color='white', fontweight='bold')

# Plot 2: Prime counting function vs approximation
ax = axes[1]
ns = np.arange(2, 1000)
p_all = sieve_eratosthenes(1000)
pi_n = np.array([np.searchsorted(p_all, n, side='right') for n in ns])
ax.plot(ns, pi_n, label='π(n): actual prime count', color='steelblue', linewidth=2)
ax.plot(ns, ns / np.log(ns), '--', label='n/ln(n): approximation', color='crimson', linewidth=2)
ax.set_xlabel('n')
ax.set_ylabel('Number of primes')
ax.set_title('Prime Number Theorem: π(n) ≈ n/ln(n)')
ax.legend()

plt.tight_layout()
plt.show()
print(f"Primes up to {N}: {primes}")

# --- Implementation: Primality test and factorization ---
import math

def is_prime_trial(n):
    """
    Trial division primality test.
    Time: O(√n)

    Args:
        n: int

    Returns:
        bool
    """
    if n < 2: return False
    if n == 2: return True
    if n % 2 == 0: return False
    for i in range(3, int(math.sqrt(n)) + 1, 2):
        if n % i == 0:
            return False
    return True

def prime_factorization(n):
    """
    Return prime factorization of n as dict {prime: exponent}.

    Args:
        n: int > 1

    Returns:
        dict: {p: e} meaning n = product of p^e
    """
    factors = {}
    d = 2
    while d * d <= n:
        while n % d == 0:
            factors[d] = factors.get(d, 0) + 1
            n //= d
        d += 1
    if n > 1:
        factors[n] = factors.get(n, 0) + 1
    return factors

def euclid_proof_infinitely_many_primes(n_steps=8):
    """
    Euclid's proof: given any finite set of primes, we can always find another.
    Construct: m = (product of known primes) + 1
    m is either prime or has a prime factor not in our list.
    """
    primes_found = [2]
    for _ in range(n_steps):
        m = 1
        for p in primes_found:
            m *= p
        m += 1  # m = p1*p2*...*pk + 1
        # Find smallest prime factor of m
        factors = prime_factorization(m)
        new_prime = min(factors.keys())
        if new_prime not in primes_found:
            primes_found.append(new_prime)
        else:
            primes_found.append(min(prime_factorization(m).keys()))
    return primes_found

print("Prime factorizations:")
for n in [12, 360, 1024, 1000000]:
    f = prime_factorization(n)
    expr = ' × '.join(f'{p}^{e}' if e > 1 else str(p) for p, e in f.items())
    print(f"  {n} = {expr}")

print("\nEuclid's construction (infinitely many primes):")
print(euclid_proof_infinitely_many_primes())

# --- Experiment 1: Twin primes and Goldbach's conjecture ---
# Hypothesis: twin primes (p, p+2) seem to be infinite but this is unproved
# Try changing: N

import numpy as np
N = 1000  # <-- modify this

def sieve_eratosthenes(n):
    is_prime = np.ones(n + 1, dtype=bool)
    is_prime[0] = is_prime[1] = False
    for i in range(2, int(n**0.5) + 1):
        if is_prime[i]:
            is_prime[i*i::i] = False
    return np.where(is_prime)[0]

primes = set(sieve_eratosthenes(N))
twin_primes = [(p, p+2) for p in range(3, N-1) if p in primes and (p+2) in primes]
print(f"Twin prime pairs up to {N}: {len(twin_primes)}")
print(f"First 10: {twin_primes[:10]}")

# Goldbach: every even integer > 2 = sum of two primes
violations = []
for n in range(4, min(N, 100), 2):
    found = any((n - p) in primes for p in primes if p <= n)
    if not found:
        violations.append(n)
print(f"\nGoldbach violations in [4, 100]: {violations} (expected: none)")

# --- Mini Project: Prime density and the Riemann Hypothesis ---
# Problem: The prime counting function π(n) has been the subject of the deepest
#          unsolved problem in mathematics.
# Task: Visualize π(n), its approximation n/ln(n), and Li(n) = ∫2^n dt/ln(t).

import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-whitegrid')

def sieve_eratosthenes(n):
    is_prime = np.ones(n + 1, dtype=bool)
    is_prime[0] = is_prime[1] = False
    for i in range(2, int(n**0.5) + 1):
        if is_prime[i]: is_prime[i*i::i] = False
    return np.where(is_prime)[0]

def li_approximation(n_vals):
    """Logarithmic integral Li(x) approximated numerically."""
    from scipy import integrate
    # TODO: compute Li(x) = integral from 2 to x of dt/ln(t)
    # Use a simple numerical integral if scipy unavailable:
    results = []
    for x in n_vals:
        t = np.linspace(2.01, max(x, 2.02), 1000)
        results.append(np.trapz(1/np.log(t), t))
    return np.array(results)

N = 5000
primes = sieve_eratosthenes(N)
ns = np.arange(2, N)
pi_n = np.searchsorted(primes, ns, side='right')

fig, axes = plt.subplots(1, 2, figsize=(13, 5))

ax = axes[0]
ax.plot(ns, pi_n, label='π(n)', color='steelblue', linewidth=1.5)
ax.plot(ns, ns / np.log(ns), '--', label='n/ln(n)', color='crimson', alpha=0.7)
ax.set_xlabel('n')
ax.set_ylabel('π(n)')
ax.set_title('Prime Counting Function')
ax.legend()

ax = axes[1]
ratio = pi_n / (ns / np.log(ns))
ax.plot(ns, ratio, color='purple', linewidth=1)
ax.axhline(1, color='gray', linestyle='--', linewidth=1, label='Perfect approximation')
ax.set_xlabel('n')
ax.set_ylabel('π(n) / (n/ln n)')
ax.set_title('Ratio π(n) / (n/ln n) → 1 as n → ∞')
ax.legend()

plt.tight_layout()
plt.show()