spb/prime-mystery-engine

1#!/usr/bin/env python32# =============================================================================3# model_c4.py — Cycle 2 PROVER: HL triple-correlation model for the C4 constant4# Author: Simon-Pierre Boucher — contact@spboucher.ai5# =============================================================================6# Heuristic model for the joint law of consecutive gaps (g1, g2) at scale7# lambda = ln x:8#     P(g1, g2)  PROPORTIONAL TO  W(g1, g2) * exp(-(g1+g2)/lambda)9# where W is the Hardy-Littlewood singular series of the triple {0, g1, g1+g2}10# restricted to primes p <= diameter (the tail p > g1+g2 is constant across11# pairs and cancels in the normalization):12#     W = prod_{p odd, p <= G} [ (1 - nu_p/p) / (1 - 1/p)^3 ] * [p=2 term const]13#     nu_p = #{0, g1, g1+g2 mod p}  (in {1,2,3}; equals 3 unless p | g1,14#            p | g2 or p | g1+g2)15# KNOWN LIMITATION (documented): no inclusion-exclusion for "no prime strictly16# inside the gaps" — this is the crude first-order model. If W == 1 (Cramer),17# rho == 0 identically; all predicted anticorrelation is singular-series18# coupling (the p | g1+g2 factor prevents factorization W = u(g1) v(g2)).19# Output: model rho(lambda) * lambda for lambdas of interest + large-lambda20# trend, printed and saved to data/model_c4.json. Deterministic.21# =============================================================================2223import json24from math import exp, log25from pathlib import Path2627import numpy as np2829DATA = Path(__file__).resolve().parent.parent / "data"303132def small_primes(n):33    s = np.ones(n + 1, dtype=bool); s[:2] = False34    for i in range(2, int(n ** 0.5) + 1):35        if s[i]:36            s[i * i:: i] = False37    return np.nonzero(s)[0]383940def model_rho(lam: float, gmax: int) -> float:41    """Pearson correlation of (g1,g2) under the HL-coupled exponential model."""42    gs = np.arange(2, gmax + 1, 2, dtype=np.int64)43    g1 = gs[:, None]44    g2 = gs[None, :]45    gsum = g1 + g246    logW = np.zeros((len(gs), len(gs)), dtype=np.float64)47    for p in small_primes(2 * gmax):48        if p == 2:49            continue  # constant factor for even gaps (nu=1), cancels50        # nu_p on the grid51        d1 = (g1 % p == 0)52        d2 = (g2 % p == 0)53        ds = (gsum % p == 0)54        nu = np.full(logW.shape, 3, dtype=np.int64)55        nu[d1 | d2 | ds] = 256        nu[d1 & d2] = 157        logW += np.log(1 - nu / p) - 3 * log(1 - 1 / p)58    P = np.exp(logW - gsum / lam)59    P /= P.sum()60    m1 = (P * g1).sum(); m2 = (P * g2).sum()61    v1 = (P * (g1 - m1) ** 2).sum(); v2 = (P * (g2 - m2) ** 2).sum()62    cov = (P * (g1 - m1) * (g2 - m2)).sum()63    return float(cov / (v1 * v2) ** 0.5)646566def main():67    out = {"model": "P(g1,g2) ~ W_HL(0,g1,g1+g2) * exp(-(g1+g2)/lambda)",68           "predictions": {}}69    print(f"{'lambda':>8} {'x=e^lam':>10} {'rho_model':>12} {'rho*lambda':>12}")70    for lam, label in [(log(1e8), "1e8"), (log(1e9), "1e9"), (log(1e10), "1e10"),71                       (log(4e10), "4e10"), (log(1e12), "1e12"),72                       (log(1e15), "1e15"), (log(1e20), "1e20"),73                       (log(1e30), "1e30"), (log(1e40), "1e40")]:74        gmax = int(14 * lam)          # truncation: weight e^-14 ~ 8e-775        r = model_rho(lam, gmax)76        out["predictions"][label] = {"lambda": round(lam, 4),77                                     "rho": round(r, 6),78                                     "rho_times_lambda": round(r * lam, 5)}79        print(f"{lam:8.3f} {label:>10} {r:12.6f} {r*lam:12.5f}")80    with open(DATA / "model_c4.json", "w") as f:81        json.dump(out, f, indent=2)82    print("written:", DATA / "model_c4.json")838485if __name__ == "__main__":86    main()87