problem_solvers/chebyshev_bias.py

"""
PROBLEM 7: Chebyshev's Bias
============================
More primes ≡ 3 (mod 4) than ≡ 1 (mod 4) below most x.
Rubinstein-Sarnak (1994) showed bias measure δ exists.
"""

import numpy as np
from typing import Dict


class ChebyshevBiasAnalyzer:
    def __init__(self):
        self.results = {}

    def sieve_primes(self, limit: int) -> np.ndarray:
        sieve = np.ones(limit + 1, dtype=bool)
        sieve[:2] = False
        for i in range(2, int(limit ** 0.5) + 1):
            if sieve[i]:
                sieve[i * i::i] = False
        return np.where(sieve)[0]

    def analyze(self, limit: int = 5_000_000) -> Dict:
        primes = self.sieve_primes(limit)

        # Separate by residue class
        mod4_1 = primes[primes % 4 == 1]
        mod4_3 = primes[primes % 4 == 3]
        mod3_1 = primes[primes % 3 == 1]
        mod3_2 = primes[primes % 3 == 2]

        # Running counts at sample points
        sample_x = np.logspace(1, np.log10(limit), 200).astype(int)
        sample_x = np.unique(sample_x[sample_x > 10])

        ratios_4 = []
        sign_changes_4 = []
        prev_sign = 1  # assume 3 > 1 initially
        for x in sample_x:
            n1 = np.sum(mod4_1 <= x)
            n3 = np.sum(mod4_3 <= x)
            ratio = n3 / max(n1, 1)
            ratios_4.append(ratio)
            sign = 1 if n3 > n1 else -1
            if sign != prev_sign:
                sign_changes_4.append(int(x))
            prev_sign = sign

        # Estimate δ = P(π(x;4,3) > π(x;4,1))
        # Approximate by proportion of sample points where 3 > 1
        delta_est = np.mean(np.array(ratios_4) > 1.0)

        self.results = {
            'limit': limit,
            'n_primes_mod4_1': len(mod4_1),
            'n_primes_mod4_3': len(mod4_3),
            'ratio_mod4_3_over_1': len(mod4_3) / max(len(mod4_1), 1),
            'delta_estimate': float(delta_est),
            'sign_changes_mod4': sign_changes_4[:10],
            'sample_x': sample_x.tolist(),
            'ratios_4': ratios_4,
        }
        return self.results

    def summary(self) -> str:
        r = self.results
        s = f"Chebyshev Bias Analysis (up to {r['limit']:,})\n{'='*50}\n"
        s += f"Primes ≡1 (mod 4): {r['n_primes_mod4_1']:,}\n"
        s += f"Primes ≡3 (mod 4): {r['n_primes_mod4_3']:,}\n"
        s += f"Ratio π(x;4,3)/π(x;4,1): {r['ratio_mod4_3_over_1']:.4f}\n"
        s += f"Estimated bias δ: {r['delta_estimate']:.4f}\n"
        s += f"Sign changes (first 10): {r['sign_changes_mod4'][:10]}\n"
        return s