problem_solvers/cramer_gaps.py

"""
PROBLEM 1: Cramér's Conjecture — Prime Gap Distribution
========================================================
BACKGROUND:
  Cramér (1936): lim sup g(p)/(ln p)² = 1
  Granville (1995): true constant is 2e^{-γ} ≈ 1.1229
  NOBODY HAS DISCRIMINATED BETWEEN THEM EMPIRICALLY.

HOW TO SOLVE:
  1. Compute prime gaps for primes up to 10^7
  2. Compute Cramér ratio g(p)/(ln p)² for ALL prime gaps
  3. Build histogram of ratio distribution
  4. Compare tail behavior to:
     - Cramér: P(ratio > λ) ~ e^{-λ}
     - Granville: P(ratio > λ) ~ e^{-2e^γ · λ} where 2e^γ ≈ 3.562
  5. Statistical test: which model fits the tail better?
"""

import numpy as np
from typing import Dict, List
from scipy import stats


class CramerGapAnalyzer:
    def __init__(self):
        self.results = {}
        self.primes = []
        self.gaps = []
        self.ratios = []

    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:
        print(f"  [Cramér] Sieving primes up to {limit:,}...")
        self.primes = self.sieve_primes(limit)
        self.gaps = np.diff(self.primes)
        log_primes = np.log(self.primes[1:])
        self.ratios = self.gaps / (log_primes ** 2)

        # Tail analysis
        thresholds = np.linspace(0.1, 2.0, 50)
        empirical_survival = []
        cramer_survival = []
        granville_survival = []

        two_e_gamma = 2.0 * np.exp(0.5772156649015329)  # ≈ 3.562

        for lam in thresholds:
            emp_surv = np.mean(self.ratios > lam)
            empirical_survival.append(emp_surv)
            cramer_survival.append(np.exp(-lam))
            granville_survival.append(np.exp(-two_e_gamma * lam))

        # Kolmogorov-Smirnov-like test on tail (λ > 0.5)
        tail_mask = thresholds > 0.5
        if np.sum(tail_mask) > 5:
            emp_tail = np.array(empirical_survival)[tail_mask]
            cra_tail = np.array(cramer_survival)[tail_mask]
            gra_tail = np.array(granville_survival)[tail_mask]
            # L2 tail error
            cra_error = np.sum((emp_tail - cra_tail) ** 2)
            gra_error = np.sum((emp_tail - gra_tail) ** 2)
        else:
            cra_error = gra_error = 0

        # Fit to Cramér model P(ratio > λ) ~ exp(-λ/c)
        valid = np.array(empirical_survival) > 0.001
        if np.sum(valid) > 5:
            log_emp = -np.log(np.array(empirical_survival)[valid])
            log_thresholds = thresholds[valid]
            cramer_fit = np.polyfit(log_thresholds, log_emp, 1)
            granville_fit = np.polyfit(log_thresholds, log_emp / two_e_gamma, 1) if two_e_gamma > 0 else [0, 0]
        else:
            cramer_fit = [1.0, 0.0]
            granville_fit = [1.0, 0.0]

        self.results = {
            'limit': limit,
            'n_primes': len(self.primes),
            'n_gaps': len(self.gaps),
            'max_ratio': float(np.max(self.ratios)),
            'mean_ratio': float(np.mean(self.ratios)),
            'std_ratio': float(np.std(self.ratios)),
            'max_gap': int(np.max(self.gaps)),
            'mean_gap': float(np.mean(self.gaps)),
            'thresholds': thresholds.tolist(),
            'empirical_survival': empirical_survival,
            'cramer_survival': cramer_survival,
            'granville_survival': granville_survival,
            'tail_error_cramer': float(cra_error),
            'tail_error_granville': float(gra_error),
            'better_model': 'Cramér' if cra_error < gra_error else 'Granville',
            'cramer_fit_slope': float(cramer_fit[0]),
            'granville_adjusted_slope': float(granville_fit[0]) if len(granville_fit) > 0 else 0,
        }
        return self.results

    def summary(self) -> str:
        r = self.results
        s = f"Cramér Gap Analysis (up to {r['limit']:,})\n{'='*50}\n"
        s += f"Primes: {r['n_primes']:,}, Gaps: {r['n_gaps']:,}\n"
        s += f"Max gap: {r['max_gap']}, Max ratio: {r['max_ratio']:.4f}\n"
        s += f"Mean gap: {r['mean_gap']:.2f}, Mean ratio: {r['mean_ratio']:.4f}\n"
        s += f"Tail fit: {r['better_model']} model is BETTER (lower L² error)\n"
        s += f"  Cramér tail L² error: {r['tail_error_cramer']:.6f}\n"
        s += f"  Granville tail L² error: {r['tail_error_granville']:.6f}\n"
        s += f"Note: Max observed ratio {r['max_ratio']:.4f} < 0.921 (record), insufficient to discriminate.\n"
        return s