spb/prime-mystery-engine
Public
Python 64.3%
TeX 35.7%
1#!/usr/bin/env python32# =============================================================================3# figures.py — paper figures from existing campaign data4# Author: Simon-Pierre Boucher — contact@spboucher.ai5# =============================================================================6# Generates (into paper/figs/):7# fig_scaling.pdf rho·lnx vs lnx: data, fit, confirmed predictions, HL model8# fig_race.pdf D(x) = N2-N4 lead changes (fine trace to 1e10 + checkpoints)9# fig_hist.pdf N(g, 1e8) with the mod-6 comb structure10# fig_certificate.pdf the 260-gap desert certificate map11# Palette (validated, dataviz method): blue #2a78d6, orange #eb6834, aqua #1baf7a12# Deterministic; the race trace is recomputed locally (scan to 1e10, ~40 s).13# =============================================================================1415import json16from math import log17from pathlib import Path1819import matplotlib20matplotlib.use("Agg")21import matplotlib.pyplot as plt22import numpy as np2324from core import primes_upto, sieve_segment2526ROOT = Path(__file__).resolve().parent.parent27FIGS = ROOT / "paper" / "figs"28FIGS.mkdir(exist_ok=True)29DATA = ROOT / "data"3031BLUE, ORANGE, AQUA = "#2a78d6", "#eb6834", "#1baf7a"32INK, MUTED = "#0b0b0b", "#52514e"3334plt.rcParams.update({35 "font.size": 8.5, "axes.labelsize": 9, "axes.titlesize": 9.5,36 "axes.spines.top": False, "axes.spines.right": False,37 "axes.linewidth": 0.7, "xtick.major.width": 0.7, "ytick.major.width": 0.7,38 "grid.linewidth": 0.5, "grid.alpha": 0.25, "legend.frameon": False,39 "figure.dpi": 150, "text.color": INK, "axes.edgecolor": MUTED,40 "xtick.color": MUTED, "ytick.color": MUTED, "axes.labelcolor": INK,41})424344def load_checkpoints():45 c2 = json.load(open(DATA / "cycle2_c4_4e10_M3U96a.json"))["checkpoints"]46 c4 = json.load(open(DATA / "cycle4_c4_1e13.json"))["checkpoints"]47 pts = [(float(x), v["rho1_lnx"], v["rho2_lnx"], v["D_N2_minus_N4"])48 for x, v in c2.items() if 1e8 <= float(x) < 1e10]49 pts += [(float(x), v["rho1_lnx"], v["rho2_lnx"], v["D_N2_minus_N4"])50 for x, v in c4.items()]51 return sorted(pts)525354def fit_cd(lams, vals):55 A = np.vstack([np.ones_like(lams), 1 / lams]).T56 (c, d), *_ = np.linalg.lstsq(A, vals, rcond=None)57 return c, d585960def fig_scaling(pts):61 lam = np.array([log(x) for x, *_ in pts])62 v1 = np.array([p[1] for p in pts])63 v2 = np.array([p[2] for p in pts])64 c1, d1 = fit_cd(lam, v1)65 c2_, d2_ = fit_cd(lam, v2)66 model = json.load(open(DATA / "model_c4.json"))["predictions"]67 ml = np.array([m["lambda"] for m in model.values()])68 mv = np.array([m["rho_times_lambda"] for m in model.values()])6970 fig, (ax, axb) = plt.subplots(1, 2, figsize=(6.4, 2.7),71 gridspec_kw={"width_ratios": [3, 2]})72 ll = np.linspace(18.0, 34.6, 200)73 # panel (a): data + fits + confirmed predictions74 ax.plot(ll, c1 + d1 / ll, color=BLUE, lw=1.4)75 ax.plot(ll, c2_ + d2_ / ll, color=AQUA, lw=1.4)76 ax.plot(lam, v1, "o", ms=3.4, color=BLUE, mfc="white", mew=1.1)77 ax.plot(lam, v2, "o", ms=3.4, color=AQUA, mfc="white", mew=1.1)78 # confirmed out-of-sample predictions (1e12 and 1e13 for lag1; 1e13 for lag2)79 for X, obs in ((1e12, -0.54756), (1e13, -0.54264)):80 ax.plot(log(X), obs, "D", ms=5.5, color=ORANGE, mfc="none", mew=1.4)81 ax.plot(log(1e13), -0.25247, "D", ms=5.5, color=ORANGE, mfc="none", mew=1.4)82 for X, txt in ((1e14, -0.5387), (1e15, -0.5350)):83 ax.plot(log(X), txt, "s", ms=4.5, color=MUTED, mfc="none", mew=1.0)84 ax.annotate("lag 1: $\\rho\\,\\ln x = c + d/\\ln x$", (18.4, -0.525),85 color=BLUE, fontsize=8)86 ax.annotate("lag 2", (19.8, -0.272), color=AQUA, fontsize=8)87 ax.annotate("confirmed predictions", (29.2, -0.572),88 color=ORANGE, fontsize=7.5, ha="center")89 ax.annotate("", xy=(log(1e12) + 0.1, -0.553), xytext=(28.2, -0.568),90 arrowprops=dict(arrowstyle="->", color=ORANGE, lw=0.9))91 ax.annotate("", xy=(log(1e13) - 0.1, -0.547), xytext=(30.0, -0.568),92 arrowprops=dict(arrowstyle="->", color=ORANGE, lw=0.9))93 ax.annotate("next targets", (31.6, -0.528), color=MUTED, fontsize=7.5)94 ax.set_xlabel("$\\ln x$")95 ax.set_ylabel("$\\rho(x)\\,\\ln x$")96 ax.grid(True)97 ax.set_xlim(17.8, 35.4)98 for X in (1e8, 1e10, 1e12, 1e14):99 ax.axvline(log(X), color=MUTED, lw=0.4, alpha=0.18)100 ax.text(log(X), ax.get_ylim()[1], f"$10^{{{int(round(log(X,10)))}}}$",101 fontsize=6.5, color=MUTED, ha="center", va="bottom")102 # panel (b): observed vs first-order HL model103 axb.plot(lam, v1, "o", ms=3.2, color=BLUE, mfc="white", mew=1.0, label="observed (lag 1)")104 axb.plot(ll, c1 + d1 / ll, color=BLUE, lw=1.2)105 axb.plot(ml, mv, "--", color=ORANGE, lw=1.4, label="first-order HL model")106 axb.axhline(0, color=MUTED, lw=0.6)107 axb.set_xlim(17.8, 95)108 axb.set_ylim(-0.65, 0.02)109 axb.set_xlabel("$\\lambda = \\ln x$")110 axb.grid(True)111 axb.annotate("observed", (23, -0.52), color=BLUE, fontsize=8)112 axb.annotate("first-order HL model\n($\\approx 32\\%$ of magnitude)",113 (40, -0.115), color=ORANGE, fontsize=8)114 fig.tight_layout()115 fig.savefig(FIGS / "fig_scaling.pdf")116 plt.close(fig)117 print("fig_scaling.pdf")118119120def race_trace(limit=10**10, per_decade=60):121 """D(x) sampled at ~per_decade log-spaced points per decade, exact."""122 targets = np.unique(np.round(np.logspace(6, np.log10(limit),123 int(per_decade * (np.log10(limit) - 6)))).astype(np.int64))124 base = primes_upto(int(limit ** 0.5) + 1)125 D = 0126 out_x, out_d = [], []127 prev = None128 ti = 0129 for lo in range(2, limit, 50_000_000):130 hi = min(lo + 50_000_000, limit)131 primes = sieve_segment(lo, hi, base)132 if prev is not None:133 primes = np.concatenate(([prev], primes))134 gaps = np.diff(primes)135 ends = primes[1:]136 delta = (gaps == 2).astype(np.int64) - (gaps == 4).astype(np.int64)137 run = D + np.cumsum(delta)138 while ti < len(targets) and targets[ti] < hi:139 j = np.searchsorted(ends, targets[ti], side="right") - 1140 out_x.append(int(targets[ti]))141 out_d.append(int(run[j]) if j >= 0 else D)142 ti += 1143 D = int(run[-1])144 prev = primes[-1]145 return np.array(out_x, dtype=float), np.array(out_d, dtype=float)146147148def fig_race(pts):149 xs, ds = race_trace()150 fig, ax = plt.subplots(figsize=(6.4, 2.5))151 ax.axhline(0, color=MUTED, lw=0.7)152 ax.plot(xs, ds, color=BLUE, lw=1.1)153 cx = [p[0] for p in pts if p[0] >= 2e10]154 cd = [p[3] for p in pts if p[0] >= 2e10]155 ax.plot(cx, cd, "D", ms=4.5, color=ORANGE, mfc="none", mew=1.2)156 ax.set_xscale("log")157 ax.set_yscale("symlog", linthresh=200)158 ax.set_xlim(1e6, 2e13)159 ax.set_ylim(-2.5e5, 6e4)160 ax.set_xlabel("$x$")161 ax.set_ylabel("$D(x) = N(2,x) - N(4,x)$")162 ax.grid(True, which="major")163 ax.annotate("first tie after $10^6$\n$x = 80{,}966{,}861$",164 xy=(8.1e7, 0), xytext=(2.5e6, -2600), fontsize=7.5, color=INK,165 arrowprops=dict(arrowstyle="->", color=MUTED, lw=0.8))166 ax.annotate("fine trace (every gap, sampled)", (2.6e6, 6000), color=BLUE, fontsize=8)167 ax.annotate("campaign checkpoints", (2.2e11, -12000), color=ORANGE, fontsize=8)168 fig.tight_layout()169 fig.savefig(FIGS / "fig_race.pdf")170 plt.close(fig)171 print("fig_race.pdf")172173174def fig_hist():175 rows = np.loadtxt(DATA / "gap_histogram_1e8.csv", delimiter=",", skiprows=1)176 g, cnt = rows[:, 0].astype(int), rows[:, 1]177 m = (g >= 2) & (g <= 150) & (g % 2 == 0)178 g, cnt = g[m], cnt[m]179 six = g % 6 == 0180 fig, ax = plt.subplots(figsize=(6.4, 2.5))181 ax.vlines(g[~six], 1, cnt[~six], color=MUTED, lw=1.6, alpha=0.75)182 ax.vlines(g[six], 1, cnt[six], color=BLUE, lw=1.9)183 ax.set_yscale("log")184 ax.set_xlabel("gap $g$")185 ax.set_ylabel("$N(g, 10^8)$")186 ax.grid(True, axis="y")187 ax.annotate("multiples of 6", (66, 2.2e5), color=BLUE, fontsize=8.5)188 ax.annotate("other even gaps", (80, 1.1e4), color=MUTED, fontsize=8.5)189 ax.annotate("jumping champion $g = 6$", xy=(6, 8.8e5), xytext=(16, 1.5e6),190 fontsize=7.5, color=INK,191 arrowprops=dict(arrowstyle="->", color=MUTED, lw=0.8))192 ax.set_xlim(0, 152)193 ax.set_ylim(1, 4e6)194 fig.tight_layout()195 fig.savefig(FIGS / "fig_hist.pdf")196 plt.close(fig)197 print("fig_hist.pdf")198199200def fig_certificate():201 cert = json.load(open(ROOT / "certs" / "desert_certificate.json"))202 interior = cert["interior_certificates"]203 cov_x, cov_y, tf_x, tf_y, mr_x = [], [], [], [], []204 for i_str, c in interior.items():205 i = int(i_str)206 if c["type"] == "covering_factor":207 cov_x.append(i); cov_y.append(c["q"])208 elif c["type"] == "trial_factor":209 tf_x.append(i); tf_y.append(c["q"])210 else:211 mr_x.append(i)212 fig, ax = plt.subplots(figsize=(6.4, 2.6))213 ax.plot(cov_x, cov_y, "o", ms=2.6, color=BLUE, mew=0)214 ax.plot(tf_x, tf_y, "D", ms=4.5, color=ORANGE, mfc="none", mew=1.2)215 top = max(tf_y) * 3216 ax.plot(mr_x, [top] * len(mr_x), "s", ms=5, color=AQUA, mfc="none", mew=1.3)217 ax.set_yscale("log")218 ax.set_xlabel("position $i$ in the desert ($N + i$, $1 \\leq i \\leq 259$)")219 ax.set_ylabel("certifying prime factor $q_i$")220 ax.grid(True, axis="y")221 ax.annotate("covering system (proven for every shift $t$)", (4, 230),222 color=BLUE, fontsize=8)223 ax.annotate("holes: explicit trial factor", (138, 2.6e4),224 color=ORANGE, fontsize=8)225 ax.annotate("hole: strong MR witness", (mr_x[0] + 8, top * 0.75),226 color=AQUA, fontsize=8)227 ax.set_xlim(0, 260)228 fig.tight_layout()229 fig.savefig(FIGS / "fig_certificate.pdf")230 plt.close(fig)231 print("fig_certificate.pdf")232233234if __name__ == "__main__":235 pts = load_checkpoints()236 fig_scaling(pts)237 fig_hist()238 fig_certificate()239 fig_race(pts) # last: recomputes a 1e10 scan (~40 s)240