spb/prime-mystery-engine

1#!/usr/bin/env python32# =============================================================================3# cycle2_c4.py — Cycle 2 ADVERSARY: C4 (scaled gap anticorrelation) to 10^10+4# Author: Simon-Pierre Boucher — contact@spboucher.ai5# =============================================================================6# Attacks C4: rho(x) < 0 and rho(x)*ln x -> c in [-0.60, -0.50] ?7# Single deterministic segmented scan of [2, LIMIT). Collects:8#   - lag-1 Pearson rho(x) and lag-2 Pearson rho2(x) of consecutive gaps,9#     at checkpoints {1,2,4}x10^k (streaming sums; junction gaps handled)10#   - C1/C3 continuation: jumping champion, G6(x) at checkpoints11#   - C2 race D = N2 - N4 at checkpoints (lead-change census continuation)12#   - C5 continuation: maximal gaps (validated against published table)13# Output: data/cycle2_c4_<L>.json  — bit-for-bit reproducible on any machine.14# Usage: python3 cycle2_c4.py 1e10 [hosttag]15# =============================================================================1617import json18import platform19import sys20import time21from math import log22from pathlib import Path2324import numpy as np2526from core import primes_upto, sieve_segment2728DATA = Path(__file__).resolve().parent.parent / "data"29DATA.mkdir(exist_ok=True)303132def pearson(S):33    n = S["n"]34    if n < 2:35        return 0.036    cov = S["xy"] / n - (S["x"] / n) * (S["y"] / n)37    vx = S["xx"] / n - (S["x"] / n) ** 238    vy = S["yy"] / n - (S["y"] / n) ** 239    return cov / (vx * vy) ** 0.5404142def acc(S, a, b):43    S["n"] += len(a)44    S["x"] += a.sum(); S["y"] += b.sum()45    S["xx"] += (a * a).sum(); S["yy"] += (b * b).sum()46    S["xy"] += (a * b).sum()474849def main():50    limit = int(float(sys.argv[1])) if len(sys.argv) > 1 else 10**1051    tag = sys.argv[2] if len(sys.argv) > 2 else platform.node()52    seg_size = 50_000_00053    cps = sorted(c for k in range(6, 12) for c in (10**k, 2 * 10**k, 4 * 10**k)54                 if c <= limit)55    bounds = sorted({limit, *cps, *range(2, limit, seg_size)} - {2})56    base = primes_upto(int(limit ** 0.5) + 1)5758    hist = np.zeros(6000, dtype=np.int64)59    S1 = dict(n=0, x=0.0, y=0.0, xx=0.0, yy=0.0, xy=0.0)   # lag-160    S2 = dict(n=0, x=0.0, y=0.0, xx=0.0, yy=0.0, xy=0.0)   # lag-261    D = 062    maximal = []63    best = 064    out_cps = {}65    prev = None66    tail = []            # last up-to-2 gaps before current segment67    t0 = time.time()68    lo = 269    cps_left = list(cps)70    for hi in bounds:71        primes = sieve_segment(lo, hi, base)72        if len(primes):73            if prev is not None:74                primes = np.concatenate(([prev], primes))75            if len(primes) >= 2:76                gaps = np.diff(primes)77                hist += np.bincount(gaps, minlength=len(hist))[: len(hist)]78                D += int((gaps == 2).sum() - (gaps == 4).sum())79                gl = np.concatenate((np.array(tail, dtype=gaps.dtype), gaps))80                k = len(tail)81                if k >= 1:82                    acc(S1, gl[k - 1:-1].astype(np.float64),83                        gl[k:].astype(np.float64))84                else:  # very first segment: all internal lag-1 pairs85                    acc(S1, gl[:-1].astype(np.float64), gl[1:].astype(np.float64))86                if k >= 2:87                    acc(S2, gl[k - 2:-2].astype(np.float64),88                        gl[k:].astype(np.float64))89                elif len(gl) > 2:  # first segments: internal lag-2 pairs90                    acc(S2, gl[:-2].astype(np.float64), gl[2:].astype(np.float64))91                if int(gaps.max()) > best:92                    starts = primes[:-1]93                    for i in np.nonzero(gaps > best)[0]:94                        g, p = int(gaps[i]), int(starts[i])95                        if g > best:96                            best = g97                            maximal.append((g, p))98                tail = [int(v) for v in gaps[-2:]]99            prev = primes[-1]100        while cps_left and hi >= cps_left[0]:101            c = cps_left.pop(0)102            nz = np.nonzero(hist)[0]103            champ = int(nz[np.argmax(hist[nz])])104            g6 = 0105            for g in range(6, int(nz.max()) - 2, 6):106                if hist[g] > hist[g - 2] and hist[g] > hist[g + 2]:107                    g6 = g108                else:109                    break110            r1, r2 = pearson(S1), pearson(S2)111            out_cps[str(c)] = {112                "rho1": round(r1, 7), "rho1_lnx": round(r1 * log(c), 5),113                "rho2": round(r2, 7), "rho2_lnx": round(r2 * log(c), 5),114                "champion": champ, "G6": g6, "D_N2_minus_N4": D,115                "n_gaps": int(S1["n"] + 1),116            }117            print(f"  cp {c:.0e}: rho1*lnx={r1*log(c):+.5f} rho2*lnx={r2*log(c):+.5f} "118                  f"champ={champ} G6={g6} D={D} [{time.time()-t0:.0f}s]", flush=True)119        lo = hi120121    out = {122        "limit": limit,123        "host": tag,124        "scan_seconds": round(time.time() - t0, 1),125        "checkpoints": out_cps,126        "maximal_gaps": [{"gap": g, "after": p, "csg": round(g / log(p) ** 2, 5)}127                         for g, p in maximal],128        "method": "deterministic segmented sieve; streaming Pearson sums; "129                  "gaps attributed to END prime; no randomness",130    }131    fname = DATA / f"cycle2_c4_{limit:.0e}_{tag}.json".replace("+", "").replace("e0", "e")132    with open(fname, "w") as f:133        json.dump(out, f, indent=2)134    print("written:", fname)135136137if __name__ == "__main__":138    main()139