spb/prime-mystery-engine

1#!/usr/bin/env python32# =============================================================================3# worker_c4.py — Cycle 3: one distributed chunk of the C4' scan4# Author: Simon-Pierre Boucher — contact@spboucher.ai5# =============================================================================6# Scans gaps whose END prime lies in [lo, hi) and emits MERGEABLE partial sums:7#   S1/S2 (lag-1/lag-2 Pearson streaming sums), gap histogram, D = N2-N4,8#   chunk-local record gaps. A 5000-wide sieve buffer below lo recovers the9#   two gaps preceding the first counted one (max gap < 1e12 is ~540, and any10#   540-window below 1e12 contains a prime, so the buffer always suffices).11# Deterministic. Usage: python3 worker_c4.py <lo> <hi> <out.json>12# =============================================================================1314import json15import platform16import sys17import time1819import numpy as np2021from core import primes_upto, sieve_segment2223import os24BUFFER = 500025SEG = int(os.environ.get("C4_SEG", 50_000_000))262728def run(lo, hi, out_path, sieve_fn=None, tag=""):29    sieve = sieve_fn or sieve_segment30    t0 = time.time()31    base = primes_upto(int(hi ** 0.5) + 1)32    start = 2 if lo <= 2 else lo - BUFFER3334    hist = np.zeros(6000, dtype=np.int64)35    S1 = dict(n=0, x=0.0, y=0.0, xx=0.0, yy=0.0, xy=0.0)36    S2 = dict(n=0, x=0.0, y=0.0, xx=0.0, yy=0.0, xy=0.0)37    D = 038    records = []39    best = 040    n_gaps = 041    prev = None42    tail = []  # up to 2 gaps immediately preceding current segment's first gap43    s = start44    while s < hi:45        e = min(s + SEG, hi)46        primes = sieve(s, e, base)47        s = e48        if len(primes) == 0:49            continue50        if prev is not None:51            primes = np.concatenate(([prev], primes))52        if len(primes) >= 2:53            gaps = np.diff(primes)54            ends = primes[1:]55            starts = primes[:-1]56            inr = ends >= lo  # gaps counted by this worker (end prime >= lo)57            gin = gaps[inr]58            n_gaps += len(gin)59            hist += np.bincount(gin, minlength=len(hist))[: len(hist)]60            D += int((gin == 2).sum() - (gin == 4).sum())61            # pair sums: pair (g_{n-1}, g_n) / (g_{n-2}, g_n) counted at g_n62            gl = np.concatenate((np.array(tail, dtype=gaps.dtype), gaps))63            k = len(tail)64            idx = np.nonzero(inr)[0] + k          # positions of counted gaps in gl65            i1 = idx[idx >= 1]66            a = gl[i1 - 1].astype(np.float64); b = gl[i1].astype(np.float64)67            S1["n"] += len(a)68            S1["x"] += float(a.sum()); S1["y"] += float(b.sum())69            S1["xx"] += float((a * a).sum()); S1["yy"] += float((b * b).sum())70            S1["xy"] += float((a * b).sum())71            i2 = idx[idx >= 2]72            a2 = gl[i2 - 2].astype(np.float64); b2 = gl[i2].astype(np.float64)73            S2["n"] += len(a2)74            S2["x"] += float(a2.sum()); S2["y"] += float(b2.sum())75            S2["xx"] += float((a2 * a2).sum()); S2["yy"] += float((b2 * b2).sum())76            S2["xy"] += float((a2 * b2).sum())77            # chunk-local records among counted gaps78            if len(gin) and int(gin.max()) > best:79                sin = starts[inr]80                for i in np.nonzero(gin > best)[0]:81                    g, p = int(gin[i]), int(sin[i])82                    if g > best:83                        best = g84                        records.append([g, p])85            tail = [int(v) for v in gaps[-2:]]86        prev = primes[-1]8788    nz = np.nonzero(hist)[0]89    out = {90        "lo": lo, "hi": hi, "n_gaps": int(n_gaps), "D": D,91        "S1": S1, "S2": S2,92        "hist": {int(g): int(hist[g]) for g in nz},93        "records": records,94        "host": platform.node() + tag,95        "seconds": round(time.time() - t0, 1),96    }97    with open(out_path, "w") as f:98        json.dump(out, f)99    print(f"done [{lo},{hi}) on {out['host']} in {out['seconds']}s", flush=True)100101102def main():103    lo, hi = int(float(sys.argv[1])), int(float(sys.argv[2]))104    out_path = sys.argv[3]105    from pathlib import Path as _P106    if _P(out_path).exists():107        print(f"skip existing {out_path}", flush=True)108        return109    if len(sys.argv) > 4 and sys.argv[4] == "gpu":110        from gpu_sieve import sieve_segment_gpu111        run(lo, hi, out_path, sieve_fn=sieve_segment_gpu, tag="/gpu")112    else:113        run(lo, hi, out_path)114115116if __name__ == "__main__":117    main()118