spb/prime-mystery-engine

1#!/usr/bin/env python32# =============================================================================3# analyze_gaps.py — Cycle 1 EXPLORER/ADVERSARY: gap statistics vs Hardy-Littlewood4# Author: Simon-Pierre Boucher — contact@spboucher.ai5# =============================================================================6# One segmented scan up to LIMIT with checkpoints. Produces (in /data):7#   gapstats_<L>.npz : per-checkpoint gap histograms, consecutive-gap pair8#                      matrix (g_n, g_{n+1}) for gaps < 400, correlation sums9#   gapstats_<L>.json: HL-normalized ratios, champion tables, correlations10# HL model used (documented):11#   N(g, x) ~ C(g) * INT_2^x exp(-g/ln t) / ln^2 t dt,12#   C(g) = 2*C2 * prod_{p | g, p odd} (p-1)/(p-2),  C2 = 0.6601618158468696.13# Deterministic; no randomness.14# Usage: python3 analyze_gaps.py 1e8 [checkpoints default 1e6,1e7,1e8,1e9,4e9 <= limit]15# =============================================================================1617import json18import sys19import time20from math import exp, log21from pathlib import Path2223import numpy as np2425from core import primes_upto, sieve_segment2627DATA = Path(__file__).resolve().parent.parent / "data"28DATA.mkdir(exist_ok=True)2930C2 = 0.660161815846869631GMAX = 400  # pair matrix covers gaps < GMAX323334def singular(g: int) -> float:35    """C(g) = 2*C2 * prod_{p|g, p>2} (p-1)/(p-2)."""36    c = 2 * C237    d, m = 3, g38    while d * d <= m:39        if m % d == 0:40            if d > 2:41                c *= (d - 1) / (d - 2)42            while m % d == 0:43                m //= d44        d += 2 if d > 2 else 145    if m > 2:46        c *= (m - 1) / (m - 2)47    return c484950def hl_integral(g: int, x: int) -> float:51    """INT_2^x exp(-g/ln t)/ln^2 t dt, trapezoid on a log grid."""52    ts = np.exp(np.linspace(log(3), log(x), 4000))53    ys = np.exp(-g / np.log(ts)) / np.log(ts) ** 254    return float(np.trapezoid(ys, ts))555657def scan(limit: int, checkpoints, segment_size: int = 50_000_000):58    base = primes_upto(int(limit ** 0.5) + 1)59    hist = np.zeros(4000, dtype=np.int64)60    pair = np.zeros((GMAX, GMAX), dtype=np.int64)61    # streaming sums for Pearson correlation of (g_n, g_{n+1})62    S = dict(n=0, x=0.0, y=0.0, xx=0.0, yy=0.0, xy=0.0)63    cp_hists, cp_corr = {}, {}64    cps = sorted(checkpoints)65    # segment boundaries aligned on checkpoints so each checkpoint histogram66    # covers exactly [2, checkpoint)67    bounds = sorted({limit, *cps, *range(2, limit, segment_size)} - {2})68    prev = None       # last prime seen69    prev_gap = None   # gap ending at prev70    t0 = time.time()71    lo = 272    for hi in bounds:73        primes = sieve_segment(lo, hi, base)74        if len(primes) == 0:75            continue76        if prev is not None:77            primes = np.concatenate(([prev], primes))78        if len(primes) >= 2:79            gaps = np.diff(primes)80            hist += np.bincount(gaps, minlength=len(hist))[: len(hist)]81            # consecutive pairs, including the junction pair across segments82            if prev_gap is not None:83                gl = np.concatenate(([prev_gap], gaps))84            else:85                gl = gaps86            a, b = gl[:-1].astype(np.float64), gl[1:].astype(np.float64)87            am, bm = gl[:-1], gl[1:]88            mask = (am < GMAX) & (bm < GMAX)89            np.add.at(pair, (am[mask], bm[mask]), 1)90            S["n"] += len(a)91            S["x"] += a.sum(); S["y"] += b.sum()92            S["xx"] += (a * a).sum(); S["yy"] += (b * b).sum()93            S["xy"] += (a * b).sum()94            prev_gap = int(gaps[-1])95        prev = primes[-1]96        while cps and hi >= cps[0]:97            c = cps.pop(0)98            cp_hists[c] = hist.copy()99            cp_corr[c] = dict(S)100            print(f"  checkpoint {c:.1e} reached at {time.time()-t0:.1f}s")101        lo = hi102    return hist, pair, cp_hists, cp_corr, time.time() - t0103104105def pearson(S):106    n = S["n"]107    cov = S["xy"] / n - (S["x"] / n) * (S["y"] / n)108    vx = S["xx"] / n - (S["x"] / n) ** 2109    vy = S["yy"] / n - (S["y"] / n) ** 2110    return cov / (vx * vy) ** 0.5111112113def main():114    limit = int(float(sys.argv[1])) if len(sys.argv) > 1 else 10**8115    cps = [int(c) for c in (10**6, 10**7, 10**8, 10**9, 4 * 10**9) if c <= limit]116    tag = f"{limit:.0e}".replace("+0", "").replace("+", "")117    hist, pair, cp_hists, cp_corr, dt = scan(limit, cps)118119    np.savez_compressed(120        DATA / f"gapstats_{tag}.npz",121        hist=hist, pair=pair,122        checkpoints=np.array(sorted(cp_hists)),123        **{f"hist_{c}": h for c, h in cp_hists.items()},124    )125126    out = {"limit": limit, "scan_seconds": round(dt, 1), "checkpoints": {}}127    for c in sorted(cp_hists):128        h = cp_hists[c]129        nz = np.nonzero(h)[0]130        even = nz[nz % 2 == 0]131        # HL comparison for even gaps up to 120132        ratios = {}133        for g in [int(v) for v in even if 2 <= v <= 120]:134            pred = singular(g) * hl_integral(g, c)135            ratios[g] = round(int(h[g]) / pred, 4) if pred > 0 else None136        champ = int(nz[np.argmax(h[nz])])137        # largest G such that every multiple of 6 up to G is a strict local max138        # of the gap histogram (h[g] > h[g-2] and h[g] > h[g+2])139        mult6_localmax_upto = 0140        for g in range(6, int(even.max()) - 2, 6):141            if h[g] > h[g - 2] and h[g] > h[g + 2]:142                mult6_localmax_upto = g143            else:144                break145        out["checkpoints"][str(c)] = {146            "jumping_champion": champ,147            "N2": int(h[2]), "N4": int(h[4]), "N6": int(h[6]),148            "N2_gt_N4": bool(h[2] > h[4]),149            "mult6_local_max_up_to": int(mult6_localmax_upto),150            "pearson_consecutive_gaps": round(pearson(cp_corr[c]), 5),151            "HL_ratio_obs_over_pred": ratios,152        }153    with open(DATA / f"gapstats_{tag}.json", "w") as f:154        json.dump(out, f, indent=2)155    # console digest156    for c in sorted(cp_hists):157        o = out["checkpoints"][str(c)]158        r = o["HL_ratio_obs_over_pred"]159        rv = [v for v in r.values() if v]160        print(f"x={c:.0e}: champ={o['jumping_champion']} N2={o['N2']} N4={o['N4']} "161              f"N2>N4={o['N2_gt_N4']} 6|g-localmax-upto={o['mult6_local_max_up_to']} "162              f"rho={o['pearson_consecutive_gaps']} HLratio[min,max]=[{min(rv):.3f},{max(rv):.3f}]")163164165if __name__ == "__main__":166    main()167