spb/prime-mystery-engine

1#!/usr/bin/env python32# =============================================================================3# explore_gaps.py — Cycle 1 EXPLORER: prime-gap statistics up to a bound4# Author: Simon-Pierre Boucher — contact@spboucher.ai5# =============================================================================6# Fully vectorized (numpy) segmented scan of [2, LIMIT). Produces:7#   data/maximal_gaps_<L>.csv   : record (maximal) gaps with merit and CSG ratio8#   data/gap_histogram_<L>.csv  : count of each even gap g9#   data/gap_mod6_start_<L>.csv : histogram of (start prime mod 6, gap mod 6)10#   data/summary_<L>.json       : headline numbers (deterministic, no seed used)11# Method: segmented Sieve of Eratosthenes (core.py), np.diff per segment,12# np.bincount accumulation. No per-prime Python loop.13# Usage: python3 explore_gaps.py 1e814# =============================================================================1516import json17import sys18import time19from math import log20from pathlib import Path2122import numpy as np2324from core import primes_upto, sieve_segment2526DATA = Path(__file__).resolve().parent.parent / "data"27DATA.mkdir(exist_ok=True)282930def scan(limit: int, segment_size: int = 50_000_000):31    base = primes_upto(int(limit ** 0.5) + 1)32    hist = np.zeros(4000, dtype=np.int64)          # gap -> count33    mod6 = np.zeros((6, 6), dtype=np.int64)        # (p mod 6, gap mod 6) -> count34    maximal = []                                   # (gap, start prime) records35    best = 036    prev = None                                    # last prime of previous segment37    n_gaps = 038    t0 = time.time()39    for lo in range(2, limit, segment_size):40        hi = min(lo + segment_size, limit)41        primes = sieve_segment(lo, hi, base)42        if len(primes) == 0:43            continue44        if prev is not None:45            primes = np.concatenate(([prev], primes))46        if len(primes) >= 2:47            gaps = np.diff(primes)48            starts = primes[:-1]49            n_gaps += len(gaps)50            hist += np.bincount(gaps, minlength=len(hist))[: len(hist)]51            np.add.at(mod6, (starts % 6, gaps % 6), 1)52            # record gaps within this segment53            seg_best = int(gaps.max())54            if seg_best > best:55                idx = np.nonzero(gaps > best)[0]56                for i in idx:57                    g, p = int(gaps[i]), int(starts[i])58                    if g > best:59                        best = g60                        maximal.append((g, p))61        prev = primes[-1]62    dt = time.time() - t063    return hist, mod6, maximal, n_gaps, dt646566def main():67    limit = int(float(sys.argv[1])) if len(sys.argv) > 1 else 10**868    tag = f"{limit:.0e}".replace("+0", "").replace("+", "")69    hist, mod6, maximal, n_gaps, dt = scan(limit)7071    # maximal gaps with merit and Cramér-Shanks-Granville ratio72    with open(DATA / f"maximal_gaps_{tag}.csv", "w") as f:73        f.write("gap,start_prime,merit,csg_ratio\n")74        for g, p in maximal:75            f.write(f"{g},{p},{g/log(p):.6f},{g/log(p)**2:.6f}\n")7677    nz = np.nonzero(hist)[0]78    with open(DATA / f"gap_histogram_{tag}.csv", "w") as f:79        f.write("gap,count\n")80        for g in nz:81            f.write(f"{g},{hist[g]}\n")8283    with open(DATA / f"gap_mod6_start_{tag}.csv", "w") as f:84        f.write("p_mod6,gap_mod6,count\n")85        for a in range(6):86            for b in range(6):87                if mod6[a, b]:88                    f.write(f"{a},{b},{mod6[a, b]}\n")8990    # headline stats91    champion = int(nz[np.argmax(hist[nz])])92    g_arr = nz.astype(float)93    c_arr = hist[nz].astype(float)94    total = c_arr.sum()95    frac_mod6 = {r: float(c_arr[g_arr % 6 == r].sum() / total) for r in range(6)}96    gmax, pmax = maximal[-1]97    best_merit = max((g / log(p), g, p) for g, p in maximal if p > 100)98    best_csg = max((g / log(p) ** 2, g, p) for g, p in maximal if p > 100)99    summary = {100        "limit": limit,101        "n_gaps": int(n_gaps),102        "scan_seconds": round(dt, 1),103        "largest_gap": {"gap": gmax, "after_prime": pmax},104        "n_maximal_gaps": len(maximal),105        "jumping_champion": champion,106        "gap_fraction_mod6": frac_mod6,107        "best_merit": {"merit": round(best_merit[0], 5), "gap": best_merit[1], "after_prime": best_merit[2]},108        "best_csg": {"csg": round(best_csg[0], 5), "gap": best_csg[1], "after_prime": best_csg[2]},109        "method": "segmented sieve of Eratosthenes, numpy, deterministic (no seed)",110    }111    with open(DATA / f"summary_{tag}.json", "w") as f:112        json.dump(summary, f, indent=2)113    print(json.dumps(summary, indent=2))114    print(f"\nmaximal gaps ({len(maximal)}):")115    for g, p in maximal:116        print(f"  gap {g:4d} after {p:>15d}  merit {g/log(p):7.4f}  csg {g/log(p)**2:.4f}")117118119if __name__ == "__main__":120    main()121