AOL × RH — run004

Overview

Run 004 extends Runs 1–3 by adding formal goodness-of-fit tests against Exponential(1) on unfolded spacings and by comparing the zeros’ spacing spectrum to lattice spokes n(r)=3r²−r+c (c∈{2..7}) via FFT overlaps. Artifacts follow the same folder layout and naming convention as previous runs.

How to reproduce locally

Pick one environment below. All produce the same r004_* figures & metrics using Odlyzko’s zeros6.gz (first 100k zeros).

1) Get the zeros file

# macOS / Linux
curl -L -o zeros6.gz "https://www-users.cse.umn.edu/~odlyzko/zeta_tables/zeros6.gz"

# Windows (PowerShell)
Invoke-WebRequest -Uri "https://www-users.cse.umn.edu/~odlyzko/zeta_tables/zeros6.gz" -OutFile "zeros6.gz"

2) Choose one runtime

Python (recommended)

1. Copy the “Reproducibility — Python (Run 004)” script from this page into run004_make_all.py.
2. Create a venv and run the script:

python -m venv .venv
source .venv/bin/activate   # Windows: .venv\Scripts\activate
pip install numpy pandas matplotlib
python run004_make_all.py

PowerShell (native, no Python)

1. Copy the “Reproducibility — PowerShell (Run 004)” block into run004_make_all.ps1.
2. Execute:

.\run004_make_all.ps1 -ZerosPath zeros6.gz

PHP (CLI)

1. Copy the “Reproducibility — PHP (Run 004)” block into run004_make_all.php.
2. Run:

php run004_make_all.php zeros6.gz

C# (.NET 6 + ScottPlot)

1. Create the project and add ScottPlot:

dotnet new console -o Run004Cs
cd Run004Cs
dotnet add package ScottPlot --version 5.0.19

2. Replace Program.cs with the “Reproducibility — C# (Run 004)” code, then run:

dotnet run -- zeros6.gz

3) Outputs (all runtimes)

Tip: if you’re seeing cached images in the browser, hard-refresh or append ?v=timestamp to the image URL.

Run 004 — Summary

Notes

• KS/CvM are against Exp(1) on unfolded Δγ; extremely small p-value underflows to 0 in double-precision.
• Spectral correlation uses Pearson on log(1+magnitude) after interpolating onto zeros’ frequency grid (skip f=0).
• Best (least-negative) lattice overlap is c=7, spec-corr≈−0.097 at ≤1e6.

Methods & Parameters

• Zero set: Odlyzko’s zeros6.gz, first 100,000 ordinates.
• Spacings: nearest-neighbor Δγ; unfolded to mean 1.0 by dividing by raw mean.
• GOF tests: one-sample KS and Cramér–von Mises vs CDF F(x)=1−e−x.
• Lattice spokes: n(r)=3r²−r+c, c∈{2..7}; primes ≤ N_max=1,000,000 via sieve.
• Spectral comparison: magnitude of rFFT on zero-mean spacings; Pearson corr on log(1+|FFT|) after interpolating lattice spectrum onto zeros’ frequency grid (excluding f=0).

Results

• Unfolded spacings: mean = 1.00000, σ ≈ 0.42925, var ≈ 0.18426 (n = 99,999).
• KS vs Exp(1): D ≈ 0.291643, p ≈ 0.0.
• CvM vs Exp(1): T ≈ 3050.068.
• Lattice spokes (≤1e6):
  • c=3 — 88 primes, 87 gaps; KS D≈0.1478 (p≈0.040), CvM≈0.5829, spec-corr≈−0.303.
  • c=5 — 84 primes, 83 gaps; KS D≈0.0929 (p≈0.453), CvM≈0.1955, spec-corr≈−0.324.
  • c=7 — 103 primes, 102 gaps; KS D≈0.1590 (p≈0.010), CvM≈0.4156, spec-corr≈−0.097.

Discussion

The unfolded zeros strongly reject Poisson spacing (as expected) under KS and CvM. First-pass spectral overlaps with AOL spokes at ≤1e6 are weak/negative; the “best” spoke (c=7) is still slightly negative. These results are consistent with non-Poisson repulsion for zeros and do not provide evidence of a simple lattice-spoke spectral match at this scale.

Limitations

• Spoke samples are small (≈80–100 spacings), limiting spectral power.
• Simple periodogram correlation is sensitive to interpolation and leakage; coherence or band-specific statistics would be more stable.
• GOF vs Exp(1) is a sanity check; a direct comparison to GUE/Wigner is planned for a follow-on run.

Metrics — inline preview

KS/CvM summary (JSON)

Loading…

Lattice spokes metrics (CSV)

Reproducibility — Python (end-to-end)

Pipeline mirrors Run 002/003; only run id and filenames change to r004_*.

# Key differences for r004:
# run_base = "run004"; output filenames use r004_*.
# Additional steps: KS/CvM vs Exp(1); lattice spokes (c=2..7) primes ≤ 1e6; FFT overlap correlation.
# (You can reuse the r003 code blocks, swapping run_base and adding the new analyses.)

Reproducibility — Python (Run 004)

Generates ECDF vs Exp(1), Q–Q vs Exp(1), FFT overlap bar, FFT overlay (best spoke), and all CSV/JSON metrics.

Environment

python -m venv .venv
source .venv/bin/activate  # Windows: .venv\Scripts\activate
pip install numpy pandas matplotlib

Script

import gzip, os, json, math, numpy as np, pandas as pd
import matplotlib.pyplot as plt

run_base = "run004"
zeros_path = "zeros6.gz"

fig_dir = os.path.join(run_base, "figures")
met_dir = os.path.join(run_base, "metrics")
data_dir = os.path.join(run_base, "data")
os.makedirs(fig_dir, exist_ok=True)
os.makedirs(met_dir, exist_ok=True)
os.makedirs(data_dir, exist_ok=True)

# ---- Load zeros and compute unfolded spacings
vals = []
with gzip.open(zeros_path, "rt", encoding="utf-8", errors="ignore") as f:
    for line in f:
        line = line.strip()
        if not line: continue
        tok = line.split()[0]
        try: v = float(tok)
        except: continue
        vals.append(v)
        if len(vals) >= 120_000: break

assert len(vals) >= 100_000, f"Only {len(vals)} zeros parsed"
z100k = np.array(vals[:100_000], float)
sp_raw = np.diff(z100k)
mean_raw = float(sp_raw.mean())
sp_unf = sp_raw / mean_raw

# ---- KS/CvM vs Exp(1)
def ks_one_sample_exp1(x):
    xs = np.sort(x); n = xs.size
    F = 1.0 - np.exp(-xs)
    i = np.arange(1, n+1)
    D = max(np.max(i/n - F), np.max(F - (i-1)/n))
    t = (np.sqrt(n) + 0.12 + 0.11/np.sqrt(n)) * D
    p = 0.0
    for k in range(1, 101):
        p += (-1)**(k-1) * math.exp(-2*(k*k)*(t*t))
    p = max(0.0, min(1.0, 2*p))
    return float(D), float(p)

def cvm_against_exp1(x):
    xs = np.sort(x); n = xs.size
    F = 1.0 - np.exp(-xs)
    u = (2*np.arange(1, n+1)-1)/(2.0*n)
    T = 1.0/(12*n) + np.sum((F-u)**2)
    return float(T)

ks_D, ks_p = ks_one_sample_exp1(sp_unf)
cvm_T = cvm_against_exp1(sp_unf)

# ---- Figures: ECDF vs Exp(1) and Q-Q vs Exp(1)
xs = np.sort(sp_unf); ys = np.arange(1, xs.size+1)/xs.size
Fx = 1.0 - np.exp(-xs)
plt.figure(figsize=(7,4.2), dpi=150)
plt.step(xs, ys, where="post", label="ECDF (unfolded Δγ)")
plt.plot(xs, Fx, label="Exponential(1) CDF")
plt.xlabel("x"); plt.ylabel("CDF"); plt.title("ECDF vs Exponential(1) — Run 004")
plt.legend(); plt.tight_layout()
plt.savefig(os.path.join(fig_dir, "r004_ecdf_vs_exp.png")); plt.close()

q = np.linspace(0.001, 0.999, 2000)
exp_q = -np.log(1.0 - q)
samp_q = np.quantile(sp_unf, q)
plt.figure(figsize=(7,4.2), dpi=150)
plt.plot(exp_q, samp_q)
plt.plot(exp_q, exp_q, ls="--")
plt.xlabel("Exp(1) quantiles"); plt.ylabel("Sample quantiles (unfolded Δγ)")
plt.title("Q–Q Plot vs Exponential(1) — Run 004")
plt.tight_layout()
plt.savefig(os.path.join(fig_dir, "r004_qq_exp.png")); plt.close()

# ---- Lattice spokes primes <= 1e6: n(r)=3r^2 - r + c, c in {2..7}
Nmax = 1_000_000
is_prime = np.ones(Nmax+1, dtype=bool); is_prime[:2] = False
for p in range(2, int(Nmax**0.5)+1):
    if is_prime[p]: is_prime[p*p:Nmax+1:p] = False

def spoke_vals(c):
    out=[]; r=0
    while True:
        n = 3*r*r - r + c
        if n > Nmax: break
        if n >= 2 and is_prime[n]: out.append(n)
        r += 1
    return np.array(out, dtype=np.int64)

def safe_spacings(a):
    if a.size < 3: return np.array([], float)
    d = np.diff(a.astype(float))
    return d[(d>0) & np.isfinite(d)]

def fft_mag(a):
    a0 = a - np.mean(a); F = np.fft.rfft(a0)
    return np.abs(F), np.fft.rfftfreq(a.size, d=1.0)

# Zeros spectrum
mag_z, f_z = fft_mag(sp_unf)
pd.DataFrame({"freq":f_z, "mag":mag_z}).to_csv(os.path.join(data_dir, "r004_fft_zeros.csv"), index=False)

rows=[]; corrs=[]; c_list=range(2,8)
for c in c_list:
    vals = spoke_vals(c)
    sp = safe_spacings(vals)
    if sp.size==0:
        rows.append(dict(c=c, n_primes=int(vals.size), n_spacings=0, mean=np.nan,
                         ks_D=np.nan, ks_p=np.nan, cvm_T=np.nan, spec_corr=np.nan))
        corrs.append(np.nan); continue
    sp_unf_l = sp / sp.mean()
    D_l, p_l = ks_one_sample_exp1(sp_unf_l)
    T_l = cvm_against_exp1(sp_unf_l)
    mag_l, f_l = fft_mag(sp_unf_l)
    # interpolate onto zeros' freqs (skip 0)
    fz = f_z[1:]; mz = mag_z[1:]
    ml = np.interp(fz, f_l, mag_l, left=np.nan, right=np.nan)
    ok = np.isfinite(ml) & np.isfinite(mz)
    if ok.sum() >= 10:
        a = np.log1p(mz[ok]); b = np.log1p(ml[ok])
        spec_corr = float(np.corrcoef(a,b)[0,1]) if a.std()!=0 and b.std()!=0 else np.nan
    else:
        spec_corr = np.nan
    rows.append(dict(c=c, n_primes=int(vals.size), n_spacings=int(sp.size),
                     mean=float(sp_unf_l.mean()), ks_D=D_l, ks_p=p_l,
                     cvm_T=T_l, spec_corr=spec_corr))
    corrs.append(spec_corr)
    pd.DataFrame({"freq":f_l,"mag":mag_l}).to_csv(os.path.join(data_dir, f"r004_fft_spoke_c{c}.csv"), index=False)

df = pd.DataFrame(rows)
df.to_csv(os.path.join(met_dir, "r004_lattice_metrics.csv"), index=False)
with open(os.path.join(met_dir, "r004_ks_cvm_summary.json"), "w") as f:
    json.dump({"raw_mean_spacing_first100k":mean_raw,"ks_D_vs_exp1":ks_D,"ks_p_vs_exp1":ks_p,"cvm_T_vs_exp1":cvm_T}, f, indent=2)

# FFT overlap bar
plt.figure(figsize=(7,4.2), dpi=150)
labels=[str(c) for c in c_list]; vals=[np.nan if v is None else v for v in corrs]
plt.bar(labels, vals); plt.xlabel("Spoke parameter c")
plt.ylabel("Spectral correlation with zeros (log-magnitude)")
plt.title("FFT overlap: zeros vs lattice spokes — Run 004")
plt.tight_layout()
plt.savefig(os.path.join(fig_dir, "r004_fft_corr_bar.png")); plt.close()

# Overlay best
if np.any(np.isfinite(corrs)):
    c_best = list(c_list)[int(np.nanargmax(corrs))]
    spec_best = pd.read_csv(os.path.join(data_dir, f"r004_fft_spoke_c{c_best}.csv"))
    plt.figure(figsize=(7,4.2), dpi=150)
    plt.loglog(f_z[1:], mag_z[1:], label="Zeros Δγ (unfolded)")
    plt.loglog(spec_best["freq"].values[1:], spec_best["mag"].values[1:], label=f"Spoke c={c_best} Δn (unfolded)")
    plt.xlabel("frequency"); plt.ylabel("|FFT|")
    plt.title(f"FFT magnitude overlay — best spoke c={c_best}")
    plt.legend(); plt.tight_layout()
    plt.savefig(os.path.join(fig_dir, "r004_fft_magnitude_overlay_best.png")); plt.close()

print("Run004 complete.")

Reproducibility — PowerShell (Run 004, native)

Windows PowerShell 5+/PowerShell 7+. Uses .NET Charting for PNGs; naive DFT for spectra (dependency-free).

# run004_make_all.ps1
param(
  [string]$ZerosPath = "zeros6.gz",
  [string]$RunBase = "run004",
  [int]$Take = 100000,
  [int]$ParseCap = 120000,
  [int]$Nmax = 1000000
)
$ErrorActionPreference = "Stop"
Add-Type -AssemblyName System.IO.Compression.FileSystem
Add-Type -AssemblyName System.Windows.Forms
Add-Type -AssemblyName System.Windows.Forms.DataVisualization

# --- IO helpers
function Ensure-Dir($path){ $d=[IO.Path]::GetDirectoryName($path); if($d -and -not(Test-Path $d)){ New-Item -ItemType Directory $d | Out-Null } }
function Write-CSV($path,[double[]]$series){
  Ensure-Dir $path
  $sw = New-Object IO.StreamWriter($path,$false,[Text.Encoding]::UTF8)
  $sw.WriteLine("index,spacing")
  for($i=0;$i -lt $series.Length;$i++){ $sw.WriteLine("$i,$($series[$i].ToString('R'))") }
  $sw.Close()
}
function Save-Chart($chart,$path){ Ensure-Dir $path; $chart.SaveImage($path, 'Png') }

# --- Parse zeros (first ~120k lines)
function Read-Zeros {
  param([string]$gzPath, [int]$cap)
  $fs = [IO.File]::OpenRead($gzPath)
  $gz = New-Object IO.Compression.GZipStream($fs,[IO.Compression.CompressionMode]::Decompress)
  $sr = New-Object IO.StreamReader($gz)
  $vals = New-Object System.Collections.Generic.List[double]
  while(-not $sr.EndOfStream -and $vals.Count -lt $cap){
    $line = $sr.ReadLine().Trim()
    if($line.Length -gt 0){
      $tok = $line.Split((@(' ', "`t")), [StringSplitOptions]::RemoveEmptyEntries)[0]
      [double]$v = 0; if([double]::TryParse($tok,[ref]$v)){ $vals.Add($v) }
    }
  }
  $sr.Close(); $gz.Close(); $fs.Close()
  ,$vals.ToArray()
}

# --- math helpers
function Diff([double[]]$a){ $n=$a.Length-1; $d=New-Object double[] $n; for($i=0;$i -lt $n;$i++){ $d[$i]=$a[$i+1]-$a[$i] } ,$d }
function Unfold([double[]]$s){ $m = ($s | Measure-Object -Average).Average; if($m -eq 0 -or [double]::IsNaN($m)){ throw "Invalid mean." }
  $u = New-Object double[] $s.Length; for($i=0;$i -lt $s.Length;$i++){ $u[$i]=$s[$i]/$m }; @{ unfolded=$u; mean=$m } }
function KS-Exp1([double[]]$x){
  $xs = $x | Sort-Object; $n=$xs.Length; $F = $xs | ForEach-Object { 1.0 - [math]::Exp(-$_) }
  $Dplus = 0.0; $Dminus = 0.0
  for($i=1;$i -le $n;$i++){
    $Dplus  = [math]::Max($Dplus,  $i/$n - $F[$i-1])
    $Dminus = [math]::Max($Dminus, $F[$i-1] - ($i-1)/$n)
  }
  $D = [math]::Max($Dplus,$Dminus)
  $t = ([math]::Sqrt($n) + 0.12 + 0.11/[math]::Sqrt($n)) * $D
  $p = 0.0; for($k=1;$k -le 100;$k++){ $p += ((-1.0)**($k-1)) * [math]::Exp(-2.0 * $k*$k * $t*$t) }; $p = [math]::Max(0.0,[math]::Min(1.0,2*$p))
  return @{ D=$D; p=$p }
}
function CvM-Exp1([double[]]$x){
  $xs = $x | Sort-Object; $n=$xs.Length
  $F = $xs | ForEach-Object { 1.0 - [math]::Exp(-$_) }
  $T = 1.0/(12*$n)
  for($i=1;$i -le $n;$i++){ $u=(2*$i-1)/(2.0*$n); $d=$F[$i-1]-$u; $T += $d*$d }
  return $T
}
function DFTMag([double[]]$a){
  $N = $a.Length; $mean = ($a | Measure-Object -Average).Average
  $x = $a | ForEach-Object { $_ - $mean }
  $M = [int][math]::Floor($N/2.0)+1
  $mag = New-Object double[] $M
  for($k=0;$k -lt $M;$k++){
    $re=0.0; $im=0.0
    for($n=0;$n -lt $N;$n++){
      $ang = -2.0*[math]::PI*$k*$n/$N
      $re += $x[$n]*[math]::Cos($ang); $im += $x[$n]*[math]::Sin($ang)
    }
    $mag[$k] = [math]::Sqrt($re*$re + $im*$im)
  }
  return ,$mag
}

# --- Sieve
function New-Sieve([int]$N){
  $arr = New-Object bool[] ($N+1); for($i=0;$i -lt $arr.Length;$i++){ $arr[$i]=$true }; $arr[0]=$false; $arr[1]=$false
  for($p=2; $p*$p -le $N; $p++){ if($arr[$p]){ for($q=$p*$p; $q -le $N; $q+=$p){ $arr[$q]=$false } } }
  return ,$arr
}
function Spoke-Primes([int]$c,[bool[]]$isPrime,[int]$N){
  $vals = New-Object System.Collections.Generic.List[int]
  $r=0; while($true){
    $n = 3*$r*$r - $r + $c
    if($n -gt $N){ break }
    if($n -ge 2 -and $isPrime[$n]){ $vals.Add($n) }
    $r++
  }
  ,$vals.ToArray()
}
function Spacings([int[]]$arr){
  if($arr.Length -lt 3){ return ,@() }
  $d = New-Object double[] ($arr.Length-1)
  for($i=0;$i -lt $arr.Length-1;$i++){ $d[$i] = [double]($arr[$i+1]-$arr[$i]) }
  ,($d | Where-Object { $_ -gt 0 -and [double]::IsFinite($_) })
}

# --- prepare dirs
$figDir = Join-Path $RunBase "figures"
$metDir = Join-Path $RunBase "metrics"
$dataDir = Join-Path $RunBase "data"
New-Item -ItemType Directory $figDir -Force | Out-Null
New-Item -ItemType Directory $metDir -Force | Out-Null
New-Item -ItemType Directory $dataDir -Force | Out-Null

# --- load zeros
$vals = Read-Zeros -gzPath $ZerosPath -cap $ParseCap
if($vals.Length -lt $Take){ throw "Parsed only $($vals.Length) zeros" }
$z100k = $vals[0..($Take-1)]
$spRaw = Diff $z100k
$u = Unfold $spRaw
$spUnf = $u.unfolded; $meanRaw = $u.mean

# --- KS/CvM vs Exp(1)
$ks = KS-Exp1 $spUnf
$cvm = CvM-Exp1 $spUnf
"Raw mean: $([math]::Round($meanRaw,6)); KS D=$([math]::Round($ks.D,6)) p=$($ks.p); CvM T=$([math]::Round($cvm,3))" | Write-Host

# --- ECDF vs Exp(1)
$sorted = $spUnf | Sort-Object
$n = $sorted.Length
$y = 1..$n | ForEach-Object { $_ / $n }
$Fx = $sorted | ForEach-Object { 1.0 - [math]::Exp(-$_) }
$ch = New-Object System.Windows.Forms.DataVisualization.Charting.Chart
$ch.Width=800; $ch.Height=480
$area = New-Object System.Windows.Forms.DataVisualization.Charting.ChartArea
$ch.ChartAreas.Add($area)
$S1 = New-Object System.Windows.Forms.DataVisualization.Charting.Series; $S1.ChartType='FastLine'; $S1.BorderWidth=2
for($i=0;$i -lt $n;$i++){ [void]$S1.Points.AddXY($sorted[$i], $y[$i]) }
$S2 = New-Object System.Windows.Forms.DataVisualization.Charting.Series; $S2.ChartType='FastLine'; $S2.BorderWidth=2
for($i=0;$i -lt $n;$i++){ [void]$S2.Points.AddXY($sorted[$i], $Fx[$i]) }
$ch.Series.Add($S1); $ch.Series.Add($S2); $ch.Titles.Add("ECDF vs Exponential(1) — Run 004") | Out-Null
Save-Chart $ch (Join-Path $figDir "r004_ecdf_vs_exp.png")

# --- Q-Q vs Exp(1)
$q = 1..2000 | ForEach-Object { ($_/2001.0) }
# simple quantiles
$idx = $q | ForEach-Object { [int][math]::Round($_*$n) - 1 } | ForEach-Object { if($_ -lt 0){0} elseif($_ -ge $n){$n-1} else {$_} }
$sampleQ = @(); foreach($i in $idx){ $sampleQ += $sorted[$i] }
$expQ = $q | ForEach-Object { -[math]::Log(1.0 - $_) }
$ch2 = New-Object System.Windows.Forms.DataVisualization.Charting.Chart
$ch2.Width=800; $ch2.Height=480
$area2 = New-Object System.Windows.Forms.DataVisualization.Charting.ChartArea
$ch2.ChartAreas.Add($area2)
$S = New-Object System.Windows.Forms.DataVisualization.Charting.Series; $S.ChartType='FastLine'; $S.BorderWidth=2
for($i=0;$i -lt $sampleQ.Count;$i++){ [void]$S.Points.AddXY([double]$expQ[$i], [double]$sampleQ[$i]) }
$ch2.Series.Add($S); $ch2.Titles.Add("Q–Q vs Exponential(1) — Run 004") | Out-Null
Save-Chart $ch2 (Join-Path $figDir "r004_qq_exp.png")

# --- Spectra & lattice spokes
# zeros spectrum
$magZ = DFTMag $spUnf
# write zeros spectrum CSV (freq index as k)
$csvZ = Join-Path $dataDir "r004_fft_zeros.csv"
$sw = New-Object IO.StreamWriter($csvZ,$false,[Text.Encoding]::UTF8)
$sw.WriteLine("freq,mag"); for($k=0;$k -lt $magZ.Length;$k++){ $sw.WriteLine("$k,$($magZ[$k].ToString('R'))") }; $sw.Close()

# sieve
$isPrime = New-Sieve $Nmax

$rows = @()
$corrs = @()
$cSet = 2..7
foreach($c in $cSet){
  $valsC = Spoke-Primes -c $c -isPrime $isPrime -N $Nmax
  $sp = Spacings $valsC
  if($sp.Length -lt 3){
    $rows += [pscustomobject]@{c=$c;n_primes=$valsC.Length;n_spacings=0;mean=[double]::NaN;ks_D=[double]::NaN;ks_p=[double]::NaN;cvm_T=[double]::NaN;spec_corr=[double]::NaN}
    $corrs += [double]::NaN; continue
  }
  $unf = Unfold $sp; $spL = $unf.unfolded
  $ksL = KS-Exp1 $spL
  $cvmL = CvM-Exp1 $spL
  $magL = DFTMag $spL

  # interpolate onto zeros frequency indices (skip k=0)
  $mz = $magZ[1..($magZ.Length-1)]
  $ml = $magL[1..([Math]::Min($magL.Length-1,$mz.Length))]
  # pad/truncate
  if($ml.Length -lt $mz.Length){ $tmp = New-Object double[] $mz.Length; for($i=0;$i -lt $mz.Length;$i++){ $tmp[$i] = ($i -lt $ml.Length) ? $ml[$i] : [double]::NaN }; $ml = $tmp }

  # correlation on log1p
  $A=@(); $B=@()
  for($i=0;$i -lt $mz.Length;$i++){
    if([double]::IsNaN($ml[$i])){ continue }
    $A += [math]::Log(1.0 + $mz[$i]); $B += [math]::Log(1.0 + $ml[$i])
  }
  if($A.Count -ge 10){
    $meanA = ($A | Measure-Object -Average).Average; $meanB = ($B | Measure-Object -Average).Average
    $sa = [math]::Sqrt(($A | ForEach-Object { ($_-$meanA)*($_-$meanA) } | Measure-Object -Sum).Sum / $A.Count)
    $sb = [math]::Sqrt(($B | ForEach-Object { ($_-$meanB)*($_-$meanB) } | Measure-Object -Sum).Sum / $B.Count)
    if($sa -eq 0 -or $sb -eq 0){ $corr=[double]::NaN }
    else{
      $cov = ($A | ForEach-Object -Begin { $i=0;$s=0.0 } -Process { $s += ($_-$meanA)*($B[$i]-$meanB); $i++ } -End { $s }) / $A.Count
      $corr = $cov/($sa*$sb)
    }
  } else { $corr=[double]::NaN }

  $rows += [pscustomobject]@{c=$c;n_primes=$valsC.Length;n_spacings=$sp.Length;mean=(($spL|Measure-Object -Average).Average);ks_D=$ksL.D;ks_p=$ksL.p;cvm_T=$cvmL;spec_corr=$corr}

  # save lattice spectrum
  $csvL = Join-Path $dataDir ("r004_fft_spoke_c{0}.csv" -f $c)
  $sw = New-Object IO.StreamWriter($csvL,$false,[Text.Encoding]::UTF8)
  $sw.WriteLine("freq,mag"); for($k=0;$k -lt $magL.Length;$k++){ $sw.WriteLine("$k,$($magL[$k].ToString('R'))") }; $sw.Close()
  $corrs += $corr
}

# write lattice metrics CSV
$csvMet = Join-Path $metDir "r004_lattice_metrics.csv"
$rows | Export-Csv -Path $csvMet -NoTypeInformation

# bar figure (correlations)
$chB = New-Object System.Windows.Forms.DataVisualization.Charting.Chart
$chB.Width=800; $chB.Height=480
$areaB = New-Object System.Windows.Forms.DataVisualization.Charting.ChartArea
$chB.ChartAreas.Add($areaB)
$bar = New-Object System.Windows.Forms.DataVisualization.Charting.Series
$bar.ChartType='Column'; $bar.BorderWidth=1
for($i=0;$i -lt $cSet.Count;$i++){ [void]$bar.Points.AddXY([string]$cSet[$i], [double]$corrs[$i]) }
$chB.Series.Add($bar); $chB.Titles.Add("FFT overlap: zeros vs lattice spokes — Run 004") | Out-Null
Save-Chart $chB (Join-Path $figDir "r004_fft_corr_bar.png")

# overlay (best)
$bestIdx = -1; $best = -1e9
for($i=0;$i -lt $corrs.Count;$i++){ if([double]::IsNaN($corrs[$i])){ continue } if($corrs[$i] -gt $best){ $best=$corrs[$i]; $bestIdx=$i } }
if($bestIdx -ge 0){
  $cBest = $cSet[$bestIdx]
  $fftBest = Import-Csv (Join-Path $dataDir ("r004_fft_spoke_c{0}.csv" -f $cBest)) | ForEach-Object { [double]$_."mag" }
  $chO = New-Object System.Windows.Forms.DataVisualization.Charting.Chart
  $chO.Width=800; $chO.Height=480
  $areaO = New-Object System.Windows.Forms.DataVisualization.Charting.ChartArea
  $areaO.AxisX.IsLogarithmic=$true; $areaO.AxisY.IsLogarithmic=$true
  $chO.ChartAreas.Add($areaO)
  $sZ = New-Object System.Windows.Forms.DataVisualization.Charting.Series; $sZ.ChartType='FastLine'; $sZ.BorderWidth=2
  for($k=1;$k -lt $magZ.Length;$k++){ [void]$sZ.Points.AddXY([double]$k, [double]$magZ[$k]) }
  $sL = New-Object System.Windows.Forms.DataVisualization.Charting.Series; $sL.ChartType='FastLine'; $sL.BorderWidth=2
  for($k=1;$k -lt $fftBest.Length;$k++){ [void]$sL.Points.AddXY([double]$k, [double]$fftBest[$k]) }
  $chO.Series.Add($sZ); $chO.Series.Add($sL); $chO.Titles.Add("FFT magnitude overlay — best spoke c=$cBest") | Out-Null
  Save-Chart $chO (Join-Path $figDir "r004_fft_magnitude_overlay_best.png")
}

# KS/CvM summary JSON
$sumPath = Join-Path $metDir "r004_ks_cvm_summary.json"
@{ raw_mean_spacing_first100k = $meanRaw; ks_D_vs_exp1 = $ks.D; ks_p_vs_exp1 = $ks.p; cvm_T_vs_exp1 = $cvm } |
  ConvertTo-Json | Set-Content -Path $sumPath -Encoding UTF8

Write-Host "Run004 complete."

Reproducibility — PHP (CLI, Run 004)

CLI usage: php run004_make_all.php zeros6.gz. Outputs CSV/JSON and run004/report.html with inline SVG figures.

<?php
// run004_make_all.php
$zeros = $argv[1] ?? 'zeros6.gz';
$run = 'run004';
@mkdir("$run/metrics", 0777, true);
@mkdir("$run/figures", 0777, true); // (unused by this PHP; figures go to HTML report)
@mkdir("$run/data", 0777, true);

// --- parse zeros (~120k lines); analyze first 100k
$vals = [];
$f = gzopen($zeros, 'r'); if(!$f) die("Cannot open $zeros\n");
while(!gzeof($f) && count($vals) < 120000){
  $line = trim(gzgets($f)); if($line==='') continue;
  $tok = preg_split('/\s+/', $line)[0];
  if(is_numeric($tok)) $vals[] = (float)$tok;
}
gzclose($f);
if(count($vals) < 100000) die("Parsed only ".count($vals)." zeros\n");
$z100k = array_slice($vals, 0, 100000);
$spRaw = [];
for($i=0;$i<count($z100k)-1;$i++){ $spRaw[] = $z100k[$i+1]-$z100k[$i]; }
$meanRaw = array_sum($spRaw)/count($spRaw);
$spUnf = array_map(fn($v) => $v/$meanRaw, $spRaw);

// --- KS/CvM vs Exp(1)
function ks_exp1($x){
  sort($x); $n=count($x);
  $Dplus=0.0; $Dminus=0.0;
  for($i=1;$i<=$n;$i++){
    $F = 1.0 - exp(-$x[$i-1]);
    $Dplus = max($Dplus,  $i/$n - $F);
    $Dminus= max($Dminus, $F - ($i-1)/$n);
  }
  $D = max($Dplus,$Dminus);
  $t = (sqrt($n)+0.12+0.11/sqrt($n))*$D;
  $p=0.0; for($k=1;$k<=100;$k++){ $p += ((-1)**($k-1)) * exp(-2*$k*$k*$t*$t); }
  $p = max(0.0, min(1.0, 2.0*$p));
  return [$D,$p];
}
function cvm_exp1($x){
  sort($x); $n=count($x);
  $T = 1.0/(12.0*$n);
  for($i=1;$i<=$n;$i++){
    $F = 1.0 - exp(-$x[$i-1]);
    $u = (2*$i-1)/(2.0*$n);
    $d = $F - $u; $T += $d*$d;
  }
  return $T;
}
[$ksD,$ksP] = ks_exp1($spUnf);
$cvmT = cvm_exp1($spUnf);

// --- Sieve to 1e6
$Nmax = 1000000;
$isPrime = array_fill(0,$Nmax+1,true); $isPrime[0]=$isPrime[1]=false;
for($p=2; $p*$p<=$Nmax; $p++){ if($isPrime[$p]){ for($q=$p*$p; $q<=$Nmax; $q+=$p){ $isPrime[$q]=false; } } }

function spoke_vals($c,$isPrime,$Nmax){
  $vals=[]; $r=0;
  while(true){
    $n = 3*$r*$r - $r + $c;
    if($n > $Nmax) break;
    if($n >= 2 && $isPrime[$n]) $vals[]=$n;
    $r++;
  }
  return $vals;
}
function spacings($a){
  if(count($a)<3) return [];
  $d=[]; for($i=0;$i<count($a)-1;$i++){ $dd=$a[$i+1]-$a[$i]; if($dd>0) $d[]=(float)$dd; }
  return $d;
}
function dft_mag($arr){
  $N=count($arr); if($N<2) return [[],[]];
  $mean=array_sum($arr)/$N; $x=array_map(fn($v)=>$v-$mean,$arr);
  $M = intdiv($N,2)+1; $mag=array_fill(0,$M,0.0);
  for($k=0;$k<$M;$k++){
    $re=0.0; $im=0.0;
    for($n=0;$n<$N;$n++){
      $ang = -2.0*M_PI*$k*$n/$N;
      $re += $x[$n]*cos($ang); $im += $x[$n]*sin($ang);
    }
    $mag[$k] = sqrt($re*$re + $im*$im);
  }
  $freq = range(0,$M-1);
  return [$freq,$mag];
}

// zeros spectrum
[$fz,$mz] = dft_mag($spUnf);
$fp = fopen("$run/data/r004_fft_zeros.csv","w"); fputcsv($fp,["freq","mag"]);
for($i=0;$i<count($fz);$i++) fputcsv($fp,[$fz[$i],$mz[$i]]); fclose($fp);

// lattice spokes
$rows=[]; $corrs=[]; $cset=[2,3,4,5,6,7];
foreach($cset as $c){
  $vals = spoke_vals($c,$isPrime,$Nmax);
  $sp = spacings($vals);
  if(count($sp)<3){
    $rows[]=["c"=>$c,"n_primes"=>count($vals),"n_spacings"=>0,"mean"=>NAN,"ks_D"=>NAN,"ks_p"=>NAN,"cvm_T"=>NAN,"spec_corr"=>NAN];
    $corrs[]=NAN; continue;
  }
  $m = array_sum($sp)/count($sp);
  $spL = array_map(fn($v)=>$v/$m,$sp);
  [$D,$p] = ks_exp1($spL);
  $T = cvm_exp1($spL);
  [$fl,$ml] = dft_mag($spL);
  // align (skip 0)
  $A=[]; $B=[];
  for($i=1;$i<min(count($ml), count($mz));$i++){
    $A[] = log(1.0 + $mz[$i]); $B[] = log(1.0 + $ml[$i]);
  }
  $corr = NAN;
  if(count($A)>=10){
    $ma = array_sum($A)/count($A); $mb = array_sum($B)/count($B);
    $sa = sqrt(array_sum(array_map(fn($v)=>($v-$ma)*($v-$ma),$A))/count($A));
    $sb = sqrt(array_sum(array_map(fn($v)=>($v-$mb)*($v-$mb),$B))/count($B));
    if($sa>0 && $sb>0){
      $cov=0.0; for($i=0;$i<count($A);$i++){ $cov += ($A[$i]-$ma)*($B[$i]-$mb); } $cov /= count($A);
      $corr = $cov/($sa*$sb);
    }
  }
  $rows[]=["c"=>$c,"n_primes"=>count($vals),"n_spacings"=>count($sp),"mean"=>array_sum($spL)/count($spL),"ks_D"=>$D,"ks_p"=>$p,"cvm_T"=>$T,"spec_corr"=>$corr];

  // save spoke spectrum CSV
  $fp = fopen("$run/data/r004_fft_spoke_c$c.csv","w"); fputcsv($fp,["freq","mag"]);
  for($i=0;$i<count($fl);$i++) fputcsv($fp,[$fl[$i],$ml[$i]]); fclose($fp);
  $corrs[]=$corr;
}

// lattice metrics CSV
$fp = fopen("$run/metrics/r004_lattice_metrics.csv","w");
fputcsv($fp,["c","n_primes","n_spacings","mean","ks_D","ks_p","cvm_T","spec_corr"]);
foreach($rows as $r){ fputcsv($fp,[$r["c"],$r["n_primes"],$r["n_spacings"],$r["mean"],$r["ks_D"],$r["ks_p"],$r["cvm_T"],$r["spec_corr"]]); }
fclose($fp);

// KS/CvM summary JSON
file_put_contents("$run/metrics/r004_ks_cvm_summary.json", json_encode([
  "raw_mean_spacing_first100k"=>$meanRaw, "ks_D_vs_exp1"=>$ksD, "ks_p_vs_exp1"=>$ksP, "cvm_T_vs_exp1"=>$cvmT
], JSON_PRETTY_PRINT));

# inline SVG helpers
function ecdf_svg($sp,$title){
  sort($sp); $n=count($sp); if($n<2) return "";
  $w=800; $h=420; $pad=50;
  $xs=$sp; $ys=[]; for($i=0;$i<$n;$i++){ $ys[] = ($i+1)/$n; }
  $minX=$xs[0]; $maxX=$xs[$n-1]; if($maxX==$minX){ $maxX=$minX+1e-6; }
  $path=""; for($i=0;$i<$n;$i++){
    $sx=$pad+($w-2*$pad)*(($xs[$i]-$minX)/($maxX-$minX));
    $sy=$h-$pad-($h-2*$pad)*$ys[$i];
    $path .= ($i==0?"M":"L").round($sx,2).",".round($sy,2)." ";
  }
  return "<h3>$title</h3><svg viewBox='0 0 $w $h' width='$w' height='$h'><rect width='$w' height='$h' fill='white' stroke='#e5e7eb'/><path d='$path' fill='none' stroke='black' stroke-width='1.5'/></svg>";
}
function qq_svg($sp,$title){
  sort($sp); $n=count($sp); if($n<10) return "";
  $w=800; $h=420; $pad=50;
  $q=[]; for($i=1;$i<2000;$i++){ $q[]=$i/2000.0; }
  $sample=[]; foreach($q as $qi){ $idx=(int)round($qi*$n)-1; if($idx<0)$idx=0; if($idx>=$n)$idx=$n-1; $sample[]=$sp[$idx]; }
  $exp=array_map(fn($u)=>-log(1.0-$u), $q);
  $minX=min($exp); $maxX=max($exp); if($maxX==$minX){ $maxX=$minX+1e-6; }
  $minY=min($sample); $maxY=max($sample); if($maxY==$minY){ $maxY=$minY+1e-6; }
  $path=""; for($i=0;$i<count($q);$i++){
    $sx=$pad+($w-2*$pad)*(($exp[$i]-$minX)/($maxX-$minX));
    $sy=$h-$pad-($h-2*$pad)*(($sample[$i]-$minY)/($maxY-$minY));
    $path.=($i==0?"M":"L").round($sx,2).",".round($sy,2)." ";
  }
  return "<h3>$title</h3><svg viewBox='0 0 $w $h' width='$w' height='$h'><rect width='$w' height='$h' fill='white' stroke='#e5e7eb'/><path d='$path' fill='none' stroke='black' stroke-width='1.5'/></svg>";
}

$report = "<!doctype html><meta charset='utf-8'><title>Run004 Report</title><h1>Run 004 — Report (PHP)</h1>
<p>Raw mean: ".number_format($meanRaw,6)." | KS D: ".number_format($ksD,6)." (p=".sprintf('%.2e',$ksP).") | CvM T: ".number_format($cvmT,3)."</p>".
ecdf_svg($spUnf, "ECDF vs Exp(1)") .
qq_svg($spUnf, "Q–Q vs Exp(1)") .
"<p>Metrics: <code>$run/metrics</code>; spectra CSVs: <code>$run/data</code>.</p>";
file_put_contents("$run/report.html", $report);

echo "Done. CSV/JSON in $run/metrics; spectra in $run/data; report at $run/report.html\n";

Reproducibility — C# (.NET 6 + ScottPlot, Run 004)

Setup

dotnet new console -o Run004Cs
cd Run004Cs
dotnet add package ScottPlot --version 5.0.19

Program.cs

using System.IO.Compression;
using ScottPlot;

string zerosPath = args.Length > 0 ? args[0] : "zeros6.gz";
string runBase = "run004";
Directory.CreateDirectory(Path.Combine(runBase,"metrics"));
Directory.CreateDirectory(Path.Combine(runBase,"figures"));
Directory.CreateDirectory(Path.Combine(runBase,"data"));

// Parse zeros (~120k lines), analyze first 100k
List<double> vals = new();
using (var fs = File.OpenRead(zerosPath))
using (var gz = new GZipStream(fs, CompressionMode.Decompress))
using (var sr = new StreamReader(gz))
{
    string? line;
    while ((line = sr.ReadLine()) != null && vals.Count < 120_000)
    {
        line = line.Trim();
        if (line.Length == 0) continue;
        var tok = line.Split((char[])null, StringSplitOptions.RemoveEmptyEntries)[0];
        if (double.TryParse(tok, out double v)) vals.Add(v);
    }
}
if (vals.Count < 100_000) throw new Exception($"Parsed only {vals.Count} zeros.");
double[] z100k = vals.Take(100_000).ToArray();
double[] spRaw = z100k.Skip(1).Zip(z100k, (a,b) => a-b).Select(x=>(double)x).ToArray();
double meanRaw = spRaw.Average();
double[] spUnf = spRaw.Select(v => v/meanRaw).ToArray();

// KS/CvM vs Exp(1)
(double D, double P) KsExp1(double[] x){
    double[] xs = x.OrderBy(v=>v).ToArray();
    int n = xs.Length; double Dplus=0, Dminus=0;
    for(int i=1;i<=n;i++){
        double F = 1.0 - Math.Exp(-xs[i-1]);
        Dplus  = Math.Max(Dplus,  i/(double)n - F);
        Dminus = Math.Max(Dminus, F - (i-1)/(double)n);
    }
    double Dval = Math.Max(Dplus,Dminus);
    double t = (Math.Sqrt(n) + 0.12 + 0.11/Math.Sqrt(n)) * Dval;
    double p=0; for(int k=1;k<=100;k++) p += Math.Pow(-1,k-1) * Math.Exp(-2*k*k*t*t);
    p = Math.Max(0, Math.Min(1, 2*p));
    return (Dval, p);
}
double CvMExp1(double[] x){
    double[] xs = x.OrderBy(v=>v).ToArray();
    int n = xs.Length; double T = 1.0/(12.0*n);
    for(int i=1;i<=n;i++){
        double F = 1.0 - Math.Exp(-xs[i-1]);
        double u = (2*i-1)/(2.0*n);
        double d = F-u; T += d*d;
    }
    return T;
}
var (ksD, ksP) = KsExp1(spUnf);
double cvmT = CvMExp1(spUnf);

// ECDF vs Exp(1)
{
    double[] xs = spUnf.OrderBy(v=>v).ToArray();
    double[] ys = Enumerable.Range(1, xs.Length).Select(i=> i/(double)xs.Length).ToArray();
    double[] Fx = xs.Select(v=> 1.0 - Math.Exp(-v)).ToArray();
    var plt = new ScottPlot.Plot(800, 480);
    plt.Add.SignalXY(xs, ys);
    plt.Add.SignalXY(xs, Fx);
    plt.Title("ECDF vs Exponential(1) — Run 004");
    plt.XLabel("x"); plt.YLabel("CDF");
    plt.SavePng(Path.Combine(runBase,"figures","r004_ecdf_vs_exp.png"), 800, 480);
}

// Q–Q vs Exp(1)
{
    double[] q = Enumerable.Range(1,1999).Select(i=> i/2000.0).ToArray();
    double[] xs = spUnf.OrderBy(v=>v).ToArray();
    int n = xs.Length;
    double[] sampleQ = q.Select(u => xs[Math.Clamp((int)Math.Round(u*n)-1,0,n-1)]).ToArray();
    double[] expQ    = q.Select(u => -Math.Log(1.0-u)).ToArray();
    var plt = new ScottPlot.Plot(800, 480);
    plt.Add.SignalXY(expQ, sampleQ);
    plt.Title("Q–Q vs Exponential(1) — Run 004");
    plt.XLabel("Exp(1) quantiles"); plt.YLabel("Sample quantiles (unfolded Δγ)");
    plt.SavePng(Path.Combine(runBase,"figures","r004_qq_exp.png"), 800, 480);
}

// Naive DFT magnitude (dependency-free)
(double[] mag, double[] freq) DFTMag(double[] a){
    int N=a.Length; double mean=a.Average();
    double[] x=a.Select(v=> v-mean).ToArray();
    int M = N/2 + 1; double[] mag=new double[M]; double[] fr=new double[M];
    for(int k=0;k<M;k++){
        double re=0, im=0;
        for(int n=0;n<N;n++){
            double ang = -2.0 * Math.PI * k * n / N;
            re += x[n]*Math.Cos(ang); im += x[n]*Math.Sin(ang);
        }
        mag[k] = Math.Sqrt(re*re + im*im); fr[k]=k;
    }
    return (mag, fr);
}

// Zeros spectrum (CSV)
var (magZ, fz) = DFTMag(spUnf);
using(var sw = new StreamWriter(Path.Combine(runBase,"data","r004_fft_zeros.csv")))
{
    sw.WriteLine("freq,mag");
    for(int i=0;i<fz.Length;i++) sw.WriteLine($"{fz[i]:R},{magZ[i]:R}");
}

// Lattice spokes primes ≤ 1e6: n(r)=3r^2 - r + c, c∈{2..7}
int Nmax=1_000_000;
bool[] isPrime = Enumerable.Repeat(true, Nmax+1).ToArray(); isPrime[0]=isPrime[1]=false;
for(int p=2;p*p<=Nmax;p++) if(isPrime[p]) for(int q=p*p;q<=Nmax;q+=p) isPrime[q]=false;

int[] SpokeVals(int c){
    List<int> outv = new();
    for(int r=0;;r++){
        long n = 3L*r*r - r + c;
        if(n>Nmax) break;
        if(n>=2 && isPrime[n]) outv.Add((int)n);
    }
    return outv.ToArray();
}
double[] Spacings(int[] a){
    if(a.Length<3) return Array.Empty<double>();
    double[] d = new double[a.Length-1];
    for(int i=0;i<d.Length;i++) d[i] = a[i+1]-a[i];
    return d.Where(v=> v>0 && double.IsFinite(v)).ToArray();
}

// Compute lattice metrics and save spectra
var rows = new List<string> { "c,n_primes,n_spacings,mean,ks_D,ks_p,cvm_T,spec_corr" };
double bestCorr = double.NegativeInfinity; int bestC = -1;

foreach(int c in new int[]{2,3,4,5,6,7}){
    var valsC = SpokeVals(c);
    var sp = Spacings(valsC);
    double corr = double.NaN;
    double mean=double.NaN, D=double.NaN, Pp=double.NaN, T=double.NaN;

    if(sp.Length >= 3){
        mean = sp.Average();
        var spL = sp.Select(v=> v/mean).ToArray();
        (D,Pp) = KsExp1(spL);
        T = CvMExp1(spL);
        var (magL, fl) = DFTMag(spL);

        // log-spectrum Pearson on overlapping indices (skip k=0)
        int M = Math.Min(magZ.Length, magL.Length) - 1;
        if(M >= 10){
            double[] A = new double[M]; double[] B = new double[M];
            for(int i=0;i<M;i++){ A[i]=Math.Log(1.0+magZ[i+1]); B[i]=Math.Log(1.0+magL[i+1]); }
            double ma=A.Average(), mb=B.Average();
            double sa=Math.Sqrt(A.Select(v=>(v-ma)*(v-ma)).Average());
            double sb=Math.Sqrt(B.Select(v=>(v-mb)*(v-mb)).Average());
            if(sa>0 && sb>0){
                double cov = 0; for(int i=0;i<M;i++) cov += (A[i]-ma)*(B[i]-mb); cov/=M;
                corr = cov/(sa*sb);
                if(!double.IsNaN(corr) && corr>bestCorr){ bestCorr=corr; bestC=c; }
            }
        }

        // save lattice spectrum CSV
        using(var sw = new StreamWriter(Path.Combine(runBase,"data",$"r004_fft_spoke_c{c}.csv")))
        {
            sw.WriteLine("freq,mag");
            for(int i=0;i<fl.Length;i++) sw.WriteLine($"{i},{magL[i]:R}");
        }
    }

    rows.Add($"{c},{valsC.Length},{sp.Length},{mean:R},{D:R},{Pp:R},{T:R},{corr:R}");
}

// lattice metrics CSV
File.WriteAllLines(Path.Combine(runBase,"metrics","r004_lattice_metrics.csv"), rows);

// FFT magnitude overlay (best spoke)
if(bestC != -1){
    var lines = File.ReadAllLines(Path.Combine(runBase,"data",$"r004_fft_spoke_c{bestC}.csv")).Skip(1);
    var magBest = lines.Select(l => double.Parse(l.Split(',')[1])).ToArray();
    var plt = new ScottPlot.Plot(800, 480);
    plt.Add.SignalXY(fz.Skip(1).Select(v=>(double)v).ToArray(), magZ.Skip(1).ToArray());
    plt.Add.SignalXY(Enumerable.Range(1, magBest.Length-1).Select(i=>(double)i).ToArray(), magBest.Skip(1).ToArray());
    plt.Axes.Bottom.ScaleLog10(true); plt.Axes.Left.ScaleLog10(true);
    plt.Title($"FFT magnitude overlay — best spoke c={bestC}");
    plt.XLabel("frequency"); plt.YLabel("|FFT|");
    plt.SavePng(Path.Combine(runBase,"figures","r004_fft_magnitude_overlay_best.png"), 800, 480);
}

// KS/CvM summary JSON
File.WriteAllText(Path.Combine(runBase,"metrics","r004_ks_cvm_summary.json"),
    System.Text.Json.JsonSerializer.Serialize(new {
        raw_mean_spacing_first100k = meanRaw,
        ks_D_vs_exp1 = ksD, ks_p_vs_exp1 = ksP, cvm_T_vs_exp1 = cvmT
    }, new System.Text.Json.JsonSerializerOptions{ WriteIndented=true })
);

Console.WriteLine("Run004 complete.");

Run

dotnet run -- zeros6.gz