Current task: —

← Session 3 — Project kickoff: RH blueprint, staged workflow + archive design

#97 assistant 2026-08-17 20:51:42.127 UTC · 718→7263 tok open ⤢
thinking (12924 chars)
Task 1.3 warm-up (thread T1): reproduce m_k(1) = 1, 4/3, 2, 13/4 for k ≤ 4 and Λ₂(0;1) = 5/36 numerically; identify the limiting density. What are these moments? "m_k(1) is the k-th moment of the limiting spectral distribution of the sine-kernel Gram matrix over the sine process". Interpretation: the Gram matrix of the Gabor frame applied to the sine process — the limiting eigenvalue distribution of G̃ under the sine process (i.e., what G̃/d converges to if zeros follow sine-process statistics). m_1 = 1 (trace), m_2 = 4/3 (the HS norm — matches R(ψ₀) = 4/3! that's why Montgomery's 4/3 appears), m_3 = 2, m_4 = 13/4. What probability distribution has moments 1, 4/3, 2, 13/4? Let me think. If X ~ Uniform[0,2]: moments 2^k/(k+1): 1, 4/3, 2, 16/5 = 3.2. No: m_4 = 13/4 = 3.25 ≠ 3.2. Close but no. What about the Marčenko–Pastur? MP with ratio λ: moments 1, 1+λ, 1+3λ+λ², ... m_2 = 4/3 → λ = 1/3. m_3 = 1+3λ+λ² = 1+1+1/9 = 19/9 ≈ 2.11 ≠ 2. No. What distribution gives exactly 1, 4/3, 2, 13/4? Consider X = Y where density is... Consider the "sine-kernel Gram matrix" — the Gram matrix of vectors v_ρ under sine process. Actually the spectral distribution of the projection/restriction: the matrix G̃/d is like a compression of the sine-kernel operator. The sine kernel K(x,y) = sin(π(x−y))/(π(x−y)) is a projection (band-limiting) operator. The Gram matrix of frame vectors at point process points ≈ (2π/L)·K restricted... The limiting spectral distribution of such Gram matrices (Gabor frames at critical density over a determinantal process)... Moments 1, 4/3, 2, 13/4: let me guess X = U·2 with... uniform[0,2] moments: 1, 4/3, 2, 16/5. m_4 = 13/4 vs 16/5: differs at 4th moment. Try X = sum of two uniforms? Irwin–Hall on [0,2]: mean 1 ✔. Variance 2·(1/12)·... sum of two U[0,1]: mean 1, var 1/6, m_2 = 1+1/6 = 7/6 ≠ 4/3. No. Try beta distributions or X ~ density (2−x)/2 on [0,2] (triangular decreasing): mean = ∫x(2−x)/2 = (1/2)(2·2 − 8/3)/... = (1/2)(4−8/3)= (1/2)(4/3) = 2/3. No. Consider density ρ(x) = (1/2)(1 + something)... Let me just compute what distribution has those moments via moment problem numerics — actually simpler: maybe it's the spectral measure of the sinc-Gram: known as the "sine-kernel Gram matrix" spectral distribution... The paper says "limiting spectral distribution of the sine-kernel Gram matrix over the sine process". This might not have a closed form; the moments 1, 4/3, 2, 13/4 were computed by the authors. Our warm-up: REPRODUCE these numbers by Monte Carlo simulation: generate sine-process points (or approximate via GUE eigenvalues) in a window, build the Gram matrix of the evaluation/Gabor vectors, compute eigenvalue moments, average. Simpler deterministic approach: the Gram matrix entries are (approximately) G_kk' = (1/aL²)Σ_γ ϕ̂(γ−α_k)ϕ̂(γ−α_k') — under sine process with intensity... In the limit, by the sine-kernel correlation, E[G] ≈ integral operator with kernel ∝ sinc(α_k − α_k')... Actually the mean Gram matrix is deterministic: E[Σ_γ f(γ−α_k)f(γ−α_k')] = ∫ f(t−α_k)f(t−α_k')ρ₁(t)dt with constant intensity ρ₁ = L/2π·(1/L)·... this reduces to the Toeplitz/sinc structure. The eigenvalue distribution of the sinc Toeplitz matrix at critical density: the sinc kernel sin(πx)/(πx) has eigenvalues... the prolate spheroidal concentration! At critical density the eigenvalue distribution of the band-limiting projection restricted to an interval — the famous prolate spheroidal "step" — eigenvalues in (0,1) with a plunge region. Hmm, but moments 1, 4/3 suggest values exceeding 1 (mean 1, second moment 4/3 > 1 means distribution extends above 1). Actually maybe simpler: the paper's m_k(1) might be moments of the Gram matrix of *random* unit vectors... Let me think about what HL*(k₀) says: tr G̃ᵏ = d·m_k(1)(1+o(1)). For k=2: ‖G̃‖² = d·4/3 — matches Theorem 5.7 with ψ₀ (R=4/3). So m_2(1) = 4/3 = R(ψ₀) — consistent with sine-process (GUE) form factor. m_1 = 1 ✔ tr G̃ = d. So m_k(1) are the moments of the limiting eigenvalue distribution of G̃ when zeros are sine-process-distributed with the flat window. These are computable by Monte Carlo: sample GUE matrix eigenvalues (sine process ≈ GUE eigenvalues locally), build G̃ with flat window, eigenvalue moments. Monte Carlo plan (E-004): 1. Sample GUE(n) matrices (n ~ 500-1000), take eigenvalues in the bulk, rescale to unit mean spacing (unfold: GUE eigenvalues in [−2√n, 2√n], bulk spacing ≈ π/√n... standard: eigenvalues/(√n) in bulk have spacing π/√n·... use standard unfolding: take eigenvalues λ_i of GUE(n), use the middle portion, scale by √n/π so mean spacing = 1). 2. Choose window [0, W] with W points... Actually mimic the paper: points γ (mean spacing 1), grid α_k with spacing 2π/L = h where L ~ log-scale — but for the sine process the relevant ratio: critical density h = 2π/L. In the paper d ≈ N and grid spacing h = 2π/L ≪ spacing 1: many grid points per zero spacing (L≈5.8, h≈1.08 ≈ spacing!). Hmm at L = log(T/2π) ~ 5-6, h ≈ 1.1 while mean spacing is ~1 — so grid ≈ zero spacing. d = LT/2π ≈ T·(L/2π) ≈ N. OK. 3. Build G̃ exactly as E-003 but with GUE points instead of zeta zeros: flat window ψ₀, ϕ = mollified 1_{[−L/2,L/2]}, v_γ = (ϕ̂(γ−α_k)). With L ~ log scale... but for the sine process the "L" parameter is free: the limit L→∞ with... The paper's asymptotics: T→∞, L = log(T/2π) → ∞ slowly. The "limiting spectral distribution" is presumably the L→∞ limit of the Gabor-Gram eigenvalue distribution over the sine process at critical density h = 2π/L. To simulate: pick moderately large L (say L = 6-12), window of W ≈ d points... this needs points in [T, 2T] with T... hmm, actually we can place the "window" anywhere: take sine-process points in [0, T_w] with grid α_k = 2πk/L, k = 0..d−1, d = ⌊LT_w/2π⌋... wait but T here anchors α_0 = T. The construction is shift-invariant; set window [0, W], d = ⌊LW/2π⌋, need W large enough that L = ... no — in the paper L = log(T/2π) is tied to T. For the sine process there's no T; L is a free parameter (the window bandwidth). The limiting distribution as L, W → ∞ with d/W → 1... The paper's HL*(k₀): tr G̃ᵏ = d m_k(1)(1+o(1)) — presumably in the joint limit. So simulate: L = 8, 10, 12; W such that d ~ 1000-2000 → W = 2πd/L. Direct approach to compute the limiting moments without full simulation: the EXPECTED Gram matrix over the sine process: E[G̃] = Toeplitz with entries ∝ sinc — eigenvalue distribution of THAT (deterministic) is the prolate-spheroidal one. But tr G̃ᵏ has fluctuations; E[tr G̃ᵏ] ≠ tr(E G̃)ᵏ — the k=2 moment includes the variance (that's where 4/3 vs 1 comes from!). Indeed E‖G̃‖² = ‖EG̃‖² + Var-terms; the 4/3 = 1 + 1/3 where 1/3 comes from the sine-kernel |K(x,y)|² term (the "1−sinc²" form factor). So the moments m_k are genuinely random-matrix quantities. Monte Carlo with GUE is the way. Simplest faithful Monte Carlo: use the actual paper construction with zeta zeros themselves but... no, we want sine-process. GUE sampling: n=1000 GUE: eigenvalues, unfold (middle 60%: scale so spacing 1), take a window of W consecutive unfolded eigenvalues ~ 800-1500, build Gram with L = something... but L should grow; with W = 1200 points and L = 10: d = LW/2π ≈ 1910 > W?? Hmm d = LW/(2π): L=10 → d ≈ 1.59W — that's not the paper's d ≈ N. In the paper d = ⌊LT/2π⌋ and N(T,2T) ≈ TL/2π·(1+2log2/L−1/L...) — d ≈ N at leading order because N(T,2T) = (T/2π)(L + 2log2 − 1) + ... hmm wait: N(T,2T) = (T/2π)(log(T/2π) + 2log2) − T/2π + O(log T) = (T/2π)(L − 1 + 2log2) ≈ (T/2π)(L − 0.386). And d = LT/2π. So d/N → L/(L−0.386) ≈ 1.077 at L=5.76 → d ≈ 1.08N. OK. For the simulation, parameters: unfolded sine-process points in [0, W] (mean spacing 1, so ~W points), L free ≥ 2-ish, grid h = 2π/L, d = ⌊LW/2π⌋. With W = 600, L = 6: d = 573. G̃ = (aL²)⁻¹ VᵀV, V = ϕ̂(γ_j − α_k) matrix (600×573). Then normalized eigenvalue moments: (1/d)tr G̃ᵏ... wait "tr G̃ᵏ = d·m_k" so m_k = (1/d)tr G̃ᵏ = mean of eigenvalues^k. Compute for k=1..4, average over ~50-200 GUE samples, extrapolate L. Expected: m_1 ≈ 1, m_2 ≈ 4/3, m_3 ≈ 2, m_4 ≈ 13/4. Christoffel Λ₂(0;1): with moments m_0..m_4 (m_0=1), the Christoffel function at 0 for the moment sequence: Λ₂(0) = 1/(Σ_{i,j≤2} (M⁻¹)_{ij}·0^{i+j})... The Christoffel function Λ_m(x) = 1/Σ_{j=0}^{m} p_j(x)² (orthonormal polys w.r.t. the distribution); at x=0: Λ₂(0) = 1/(p₀(0)² + p₁(0)² + p₂(0)²). Computable from the moment matrix M = (m_{i+j})_{i,j≤2}: Λ₂(0) = (M⁻¹)_{00}·... precisely Λ_m(x) = 1/(u*M⁻¹u) with u = (1, x, x²): at x=0: Λ₂(0) = 1/(M⁻¹)_{00}. With m_0=1, m_1=1, m_2=4/3, m_3=2, m_4=13/4: M = [[1,1,4/3],[1,4/3,2],[4/3,2,13/4]]. Compute (M⁻¹)_{00}: det stuff. M⁻¹_{00} = cofactor/det = det([[4/3,2],[2,13/4]])/det(M) = (4/3·13/4 − 4)/det(M) = (13/3 − 4)/det(M) = (1/3)/det(M). det(M) = 1·(4/3·13/4−4) − 1·(13/4−8/3) + 4/3·(2−16/9) = (1/3) − (13/4−8/3) + (4/3)(2/9) = 1/3 − (39/12−32/12) + 8/27 = 1/3 − 7/12 + 8/27. LCD 108: 36/108 − 63/108 + 32/108 = 5/108. So M⁻¹_00 = (1/3)/(5/108) = 36/5. Λ₂(0) = 5/36 ✔✔ — matches the paper! And the certificate: 1 − Λ₂(0) = 31/36... hmm the paper says HL*(4) gives 13/18 = 26/36. Note 1 − Λ₂(0) = 31/36 ≠ 13/18. Hmm — maybe n₊(G̃)/d ≥ 1 − Λ_m(0) with... 1 − 5/36 = 31/36 ≈ 0.861. But paper says HL*(4) → 13/18 ≈ 0.7222 for simple-on-line. The rank bound gives rank P₁... n₊(G̃)/d ≥ 31/36 would be on-line zeros fraction; simple needs another step (the CMS bound applies to n₊; converting to N^s_0 uses the chain with charges). Whatever — 13/18 = (4 − R(4)... let me see: 13/18 = (4−4/3... no. (2 − R₂)/... With third moment the certificate improves via CMS. Don't need to fully rederive — the numeric reproduction of 5/36 from moments 1,4/3,2,13/4 is already verified by hand above! So: Λ₂(0) = 5/36 VERIFIED by hand from the moment matrix (computation above: det = 5/108, cofactor = 1/3, ratio inverse = 5/36 ✔). Also 13/18: maybe via the weight ω(m) cubic... §7.3 mentions 13/18? No, §7.2(f): HL*(4) → 13/18 for N^s_0/N. The one-sided CMS bound on n₊... Let me not over-derive; record 5/36 verified; 13/18 derivation path noted for later. What distribution has moments 1, 4/3, 2, 13/4? Try X ~ 1 − cos? Density on [0,2]: ρ(x) = (1/π)... arcsine on [0,2]: ρ = 1/(π√(x(2−x))): mean 1 ✔; m_2 = 3/2 ≠ 4/3. No. Semi-circle on [0,2] (radius 1 centered 1): mean 1, m_2 = 1 + 1/4 = 5/4 ≠ 4/3. No. Uniform[0,2]: 1, 4/3, 2, 16/5: first three match! m_4 = 13/4 = 3.25 vs 16/5 = 3.2 — so close but distinct: NOT uniform. Interesting: matches uniform to 3rd order. The distribution is a perturbation of uniform[0,2]. Monte Carlo should reproduce m_4 = 3.25 to distinguish. Simulation size: need relative accuracy ~1% on m_4. Feasible. Plan E-004: GUE Monte Carlo: - Sample GUE(400), eigenvalues, unfold using exact GUE mean spacing: for bulk eigenvalues, unfolded x_i = λ_i·√(2n)/π·... Simpler: use known GUE semicircle: density ρ(λ) = (1/2π)√(4n−λ²)... for GUE with off-diagonal variance 1, eigenvalues in [−2√n, 2√n]; bulk mean spacing at 0: 1/ρ(0) = 2π/(2√n·... ρ(0) = (1/π)√n·... Let me just take middle eigenvalues and rescale by mean local spacing: spacings of consecutive middle eigenvalues ≈ π/√n·... I'll unfold numerically: fit the cumulative. Standard: take eigenvalues λ_(0.3n..0.7n), rescale by dividing by mean spacing. - Window [0,W]: take W_unf consecutive unfolded points ≈ 600. - L ∈ {4, 6, 8}: build G̃ as in E-003 (flat window ψ₀; note: the m_k(1) presumably for ψ₀ — the paper says "sine-kernel Gram matrix" — flat window corresponds to sinc kernel ✔). - Eigenvalues of G̃, normalized moments m_k = mean(λᵏ), k=1..4. - Average over 100 samples; report mean ± std; compare to 1, 4/3, 2, 13/4. Subtlety: with L small (4-8), the mollified ϕ ≈ indicator — the "sine kernel" is the Fourier transform of the flat window: ϕ̂(ξ) ∝ sinc(Lξ/2)·L — yes sinc with parameter L/2. Also d vs N_points: in window [0,W] with W points and grid d = LW/2π: for L=6: d = 0.955W < W; fine. Performance: d ~ 573 at L=6, W=600; V is (W+d')×d — small; 100 samples × 3 L values: trivial. But wait — critical check: is the "sine-kernel Gram matrix" in the paper built with the smoothed ϕ or the sharp indicator? Sharp indicator ⟹ ϕ̂ = sinc exactly ("sine kernel" name suggests sharp). The mollification is a technicality; the limiting distribution should use the sharp cutoff: v entries = sinc((γ−α_k)L/2)... with sharp ϕ = 1_{[−L/2,L/2]}: ϕ̂(ξ) = 2sin(Lξ/2)/ξ = L·sinc(Lξ/2). Then a = ‖ϕ‖²/L = 1. G̃ = (1/L²)VᵀV with V = sinc matrix. I'll offer both but sharp is cleaner: call it the sinc-Gram matrix directly: G̃_{kk'} = (1/L²)Σ_γ sinc(L(γ−α_k)/2)sinc(L(γ−α_k')/2)·L² = Σ_γ sinc(L(γ−α_k)/2)sinc(L(γ−α_k')/2). With sinc(x) := sin(x)/x. Good—E-004 will compute both ways? Keep simple: sharp sinc (matches "sine-kernel Gram" name). Write E-004 script, run, compare. Also record the hand-verification of Λ₂(0) = 5/36 in D-010 (short derivation note) and identify density (near-uniform perturbation). Let me write the script.
Proceeding with **task 1.3 warm-up (thread T1)**. First, a quick hand result before any numerics — the Christoffel value can be verified analytically from the moment sequence: With m₀..m₄ = 1, 1, 4/3, 2, 13/4, the moment matrix M = (m_{i+j})_{i,j≤2} has det(M) = 5/108 and (M⁻¹)₀₀ = (1/3)/(5/108) = 36/5, so **Λ₂(0) = 1/(M⁻¹)₀₀ = 5/36 ✔ — matches the paper exactly.** One number already reproduced by hand. Also notable: 1, 4/3, 2 matches Uniform[0,2] moments 2^k/(k+1) to third order, but m₄ = 13/4 = 3.25 ≠ 16/5 = 3.2 — the limiting density is a *near-uniform perturbation*. Now the Monte Carlo confirmation (E-004): sample GUE spectra (bulk = sine process), build the sinc-Gram matrix exactly as in §7.2(f), measure the eigenvalue moments.
tool: Write ok
input
{
  "file_path": "D:\\Python\\Quantum\\research\\numerics\\scripts\\E-004-sine-moments.py",
  "content": "\"\"\"E-004 \u2014 Thread T1 / task 1.3 warm-up: reproduce the limiting spectral\nmoments m_k(1) = 1, 4/3, 2, 13/4 (k <= 4) of the sine-kernel Gram matrix\nover the sine process (paper \u00a77.2(f)), and the Christoffel value\nLambda_2(0) = 5/36.\n\nConstruction (flat window; \"sine-kernel Gram matrix\"):\n  points gamma: sine-process sample (GUE bulk eigenvalues, unfolded to\n  mean spacing 1) in a window of W points;\n  grid alpha_k = 2 pi k / L at critical density h = 2pi/L, k = 0..d-1,\n  d = floor(L W / 2 pi);\n  G_tilde = V^T V with V[j,k] = sinc(L(gamma_j - alpha_k)/2),\n  sinc(x) = sin(x)/x   (sharp cutoff; a = ||phi||^2/L = 1).\nMoments: m_k = (1/d) tr G_tilde^k  = mean of eigenvalues^k.\n\nHand check (already verified, D-010): from m_0..m_4 = 1,1,4/3,2,13/4 the\nmoment matrix has det = 5/108 and (M^-1)_00 = 36/5, so Lambda_2(0) = 5/36.\n\nOutputs: research/numerics/data/E-004-sine-moments.csv (+ .json)\n\"\"\"\nimport json\nimport sys\nimport time\nfrom datetime import datetime, timezone\nfrom pathlib import Path\n\nimport numpy as np\n\nREPO = Path(__file__).resolve().parents[3]\nDATA = REPO / \"research\" / \"numerics\" / \"data\"\nCSV = DATA / \"E-004-sine-moments.csv\"\nJSON = DATA / \"E-004-sine-moments.json\"\n\nRNG_SEED = 20260817\nN_GUE = 400          # GUE size per sample (bulk used)\nW = 600              # window length in unfolded points\nL_VALUES = (4.0, 6.0, 8.0)\nN_SAMPLES = 60\nTARGET = {1: 1.0, 2: 4.0 / 3.0, 3: 2.0, 4: 13.0 / 4.0}\n\n\ndef gue_bulk_points(rng) -> np.ndarray:\n    \"\"\"Unfolded GUE bulk eigenvalues (mean spacing 1), W consecutive.\"\"\"\n    A = rng.standard_normal((N_GUE, N_GUE)) + 1j * rng.standard_normal((N_GUE, N_GUE))\n    H = (A + A.conj().T) / 2.0\n    ev = np.linalg.eigvalsh(H)  # support ~ [-2 sqrt(2n), 2 sqrt(2n)] for this normalization? no:\n    # entries variance 1 (complex): semicircle radius 2*sqrt(n)\n    ev = ev / np.sqrt(N_GUE)    # now support [-2, 2], bulk density (1/pi)... unfold:\n    # semicircle density rho(x) = (1/(2 pi)) sqrt(4 - x^2); unfold via CDF\n    # take middle 60% and rescale by local mean spacing at 0: spacing(0) = 1/rho(0) = pi\n    lo, hi = int(0.2 * N_GUE), int(0.8 * N_GUE)\n    bulk = ev[lo:hi]\n    bulk = bulk * (1.0 / np.pi) * N_GUE / (0.6 * N_GUE) * len(bulk) * (np.pi / len(bulk))\n    # simpler robust unfolding: rescale so mean consecutive spacing is exactly 1\n    bulk = bulk / np.mean(np.diff(bulk))\n    mid = len(bulk) // 2\n    return bulk[mid - W // 2: mid + W // 2]\n\n\ndef sinc(x):\n    return np.sinc(x / np.pi)  # np.sinc(u) = sin(pi u)/(pi u) -> sin(x)/x\n\n\ndef main() -> int:\n    if CSV.exists() and \"--force\" not in sys.argv:\n        print(f\"cached: {CSV} exists (use --force to recompute)\")\n        return 0\n    rng = np.random.default_rng(RNG_SEED)\n    t0 = time.time()\n\n    # accumulate moments per L\n    acc = {L: {k: [] for k in (1, 2, 3, 4)} for L in L_VALUES}\n    for s in range(N_SAMPLES):\n        pts = gue_bulk_points(rng)\n        pts = pts - pts[0]\n        Wpts = pts[-1]\n        for L in L_VALUES:\n            h = 2 * np.pi / L\n            d = int(L * Wpts / (2 * np.pi))\n            alpha = h * np.arange(d)\n            Z = pts[:, None] - alpha[None, :]\n            V = sinc(L * Z / 2.0)\n            G = (V.T @ V) / L**2\n            ev = np.linalg.eigvalsh(G)\n            for k in (1, 2, 3, 4):\n                acc[L][k].append(float(np.mean(ev**k)))\n        if (s + 1) % 20 == 0:\n            print(f\"  {s+1}/{N_SAMPLES} samples ({time.time()-t0:.0f}s)\",\n                  flush=True)\n\n    rows = []\n    checks = {}\n    for L in L_VALUES:\n        row = {\"L\": L}\n        for k in (1, 2, 3, 4):\n            arr = np.array(acc[L][k])\n            m, sd = arr.mean(), arr.std()\n            row[f\"m{k}\"] = f\"{m:.5f}\"\n            row[f\"m{k}_std\"] = f\"{sd:.5f}\"\n            row[f\"m{k}_target\"] = f\"{TARGET[k]:.5f}\"\n            row[f\"m{k}_relerr\"] = f\"{abs(m / TARGET[k] - 1):.5f}\"\n        rows.append(row)\n        # generous tolerances: finite-L + finite-W + sampling effects\n        checks[f\"L={L}\"] = {\n            f\"m{k}_within_6pct\": bool(abs(np.mean(acc[L][k]) / TARGET[k] - 1) < 0.06)\n            for k in (1, 2, 3, 4)\n        }\n        print(f\"L={L}: \" + \"  \".join(\n            f\"m{k}={np.mean(acc[L][k]):.4f}+-{np.std(acc[L][k]):.4f}\"\n            f\" (target {TARGET[k]:.4f})\" for k in (1, 2, 3, 4)), flush=True)\n\n    # Christoffel from simulated moments (L=8): Lambda_2(0) = 1/(M^-1)_00\n    Lmax = L_VALUES[-1]\n    m = [1.0] + [float(np.mean(acc[Lmax][k])) for k in (1, 2, 3, 4)]\n    M = np.array([[m[i + j] for j in range(3)] for i in range(3)])\n    lam2 = 1.0 / np.linalg.inv(M)[0, 0]\n    print(f\"Christoffel Lambda_2(0) from simulated moments (L={Lmax}): \"\n          f\"{lam2:.5f}  (target 5/36 = {5/36:.5f})\")\n    checks[\"christoffel\"] = {\"within_10pct\": bool(abs(lam2 / (5 / 36) - 1) < 0.10)}\n\n    DATA.mkdir(parents=True, exist_ok=True)\n    with CSV.open(\"w\", encoding=\"utf-8\", newline=\"\") as f:\n        f.write(\",\".join(rows[0].keys()) + \"\\n\")\n        for r in rows:\n            f.write(\",\".join(str(v) for v in r.values()) + \"\\n\")\n    meta = {\n        \"experiment\": \"E-004\",\n        \"task\": \"1.3 warm-up (thread T1)\",\n        \"description\": \"Monte Carlo reproduction of sine-kernel Gram matrix \"\n                       \"spectral moments m_k(1)=1,4/3,2,13/4 and Christoffel \"\n                       \"Lambda_2(0)=5/36 (paper \u00a77.2(f))\",\n        \"script\": \"research/numerics/scripts/E-004-sine-moments.py\",\n        \"parameters\": {\"N_GUE\": N_GUE, \"W\": W, \"L_values\": list(L_VALUES),\n                       \"n_samples\": N_SAMPLES, \"seed\": RNG_SEED},\n        \"christoffel_simulated\": lam2,\n        \"checks\": checks,\n        \"elapsed_seconds\": round(time.time() - t0, 1),\n        \"timestamp_utc\": datetime.now(timezone.utc).isoformat(),\n        \"row_count\": len(rows),\n    }\n    JSON.write_text(json.dumps(meta, indent=2), encoding=\"utf-8\")\n    ok = all(all(c.values()) for c in checks.values())\n    print(f\"wrote {CSV.name}, {JSON.name} \u2014 \"\n          f\"{'ALL CHECKS PASS' if ok else 'CHECK FAILURE (inspect)'}\")\n    return 0 if ok else 1\n\n\nif __name__ == \"__main__\":\n    raise SystemExit(main())\n"
}
result 2026-08-17 20:51:42Z → 2026-08-18 01:26:18Z
File created successfully at: D:\Python\Quantum\research\numerics\scripts\E-004-sine-moments.py
raw content_json (database row)
[{"type":"thinking","thinking":"Task 1.3 warm-up (thread T1): reproduce m_k(1) = 1, 4/3, 2, 13/4 for k ≤ 4 and Λ₂(0;1) = 5/36 numerically; identify the limiting density.\n\nWhat are these moments? \"m_k(1) is the k-th moment of the limiting spectral distribution of the sine-kernel Gram matrix over the sine process\". Interpretation: the Gram matrix of the Gabor frame applied to the sine process — the limiting eigenvalue distribution of G̃ under the sine process (i.e., what G̃/d converges to if zeros follow sine-process statistics). m_1 = 1 (trace), m_2 = 4/3 (the HS norm — matches R(ψ₀) = 4/3! that's why Montgomery's 4/3 appears), m_3 = 2, m_4 = 13/4.\n\nWhat probability distribution has moments 1, 4/3, 2, 13/4? Let me think. If X ~ Uniform[0,2]: moments 2^k/(k+1): 1, 4/3, 2, 16/5 = 3.2. No: m_4 = 13/4 = 3.25 ≠ 3.2. Close but no.\n\nWhat about the Marčenko–Pastur? MP with ratio λ: moments 1, 1+λ, 1+3λ+λ², ... m_2 = 4/3 → λ = 1/3. m_3 = 1+3λ+λ² = 1+1+1/9 = 19/9 ≈ 2.11 ≠ 2. No.\n\nWhat distribution gives exactly 1, 4/3, 2, 13/4? Consider X = Y where density is... Consider the \"sine-kernel Gram matrix\" — the Gram matrix of vectors v_ρ under sine process. Actually the spectral distribution of the projection/restriction: the matrix G̃/d is like a compression of the sine-kernel operator. The sine kernel K(x,y) = sin(π(x−y))/(π(x−y)) is a projection (band-limiting) operator. The Gram matrix of frame vectors at point process points ≈ (2π/L)·K restricted... The limiting spectral distribution of such Gram matrices (Gabor frames at critical density over a determinantal process)...\n\nMoments 1, 4/3, 2, 13/4: let me guess X = U·2 with... uniform[0,2] moments: 1, 4/3, 2, 16/5. m_4 = 13/4 vs 16/5: differs at 4th moment. Try X = sum of two uniforms? Irwin–Hall on [0,2]: mean 1 ✔. Variance 2·(1/12)·... sum of two U[0,1]: mean 1, var 1/6, m_2 = 1+1/6 = 7/6 ≠ 4/3. No.\n\nTry beta distributions or X ~ density (2−x)/2 on [0,2] (triangular decreasing): mean = ∫x(2−x)/2 = (1/2)(2·2 − 8/3)/... = (1/2)(4−8/3)= (1/2)(4/3) = 2/3. No.\n\nConsider density ρ(x) = (1/2)(1 + something)... Let me just compute what distribution has those moments via moment problem numerics — actually simpler: maybe it's the spectral measure of the sinc-Gram: known as the \"sine-kernel Gram matrix\" spectral distribution... The paper says \"limiting spectral distribution of the sine-kernel Gram matrix over the sine process\". This might not have a closed form; the moments 1, 4/3, 2, 13/4 were computed by the authors. Our warm-up: REPRODUCE these numbers by Monte Carlo simulation: generate sine-process points (or approximate via GUE eigenvalues) in a window, build the Gram matrix of the evaluation/Gabor vectors, compute eigenvalue moments, average.\n\nSimpler deterministic approach: the Gram matrix entries are (approximately) G_kk' = (1/aL²)Σ_γ ϕ̂(γ−α_k)ϕ̂(γ−α_k') — under sine process with intensity... In the limit, by the sine-kernel correlation, E[G] ≈ integral operator with kernel ∝ sinc(α_k − α_k')... Actually the mean Gram matrix is deterministic: E[Σ_γ f(γ−α_k)f(γ−α_k')] = ∫ f(t−α_k)f(t−α_k')ρ₁(t)dt with constant intensity ρ₁ = L/2π·(1/L)·... this reduces to the Toeplitz/sinc structure. The eigenvalue distribution of the sinc Toeplitz matrix at critical density: the sinc kernel sin(πx)/(πx) has eigenvalues... the prolate spheroidal concentration! At critical density the eigenvalue distribution of the band-limiting projection restricted to an interval — the famous prolate spheroidal \"step\" — eigenvalues in (0,1) with a plunge region. Hmm, but moments 1, 4/3 suggest values exceeding 1 (mean 1, second moment 4/3 > 1 means distribution extends above 1).\n\nActually maybe simpler: the paper's m_k(1) might be moments of the Gram matrix of *random* unit vectors... Let me think about what HL*(k₀) says: tr G̃ᵏ = d·m_k(1)(1+o(1)). For k=2: ‖G̃‖² = d·4/3 — matches Theorem 5.7 with ψ₀ (R=4/3). So m_2(1) = 4/3 = R(ψ₀) — consistent with sine-process (GUE) form factor. m_1 = 1 ✔ tr G̃ = d. So m_k(1) are the moments of the limiting eigenvalue distribution of G̃ when zeros are sine-process-distributed with the flat window. These are computable by Monte Carlo: sample GUE matrix eigenvalues (sine process ≈ GUE eigenvalues locally), build G̃ with flat window, eigenvalue moments.\n\nMonte Carlo plan (E-004):\n1. Sample GUE(n) matrices (n ~ 500-1000), take eigenvalues in the bulk, rescale to unit mean spacing (unfold: GUE eigenvalues in [−2√n, 2√n], bulk spacing ≈ π/√n... standard: eigenvalues/(√n) in bulk have spacing π/√n·... use standard unfolding: take eigenvalues λ_i of GUE(n), use the middle portion, scale by √n/π so mean spacing = 1).\n2. Choose window [0, W] with W points... Actually mimic the paper: points γ (mean spacing 1), grid α_k with spacing 2π/L = h where L ~ log-scale — but for the sine process the relevant ratio: critical density h = 2π/L. In the paper d ≈ N and grid spacing h = 2π/L ≪ spacing 1: many grid points per zero spacing (L≈5.8, h≈1.08 ≈ spacing!). Hmm at L = log(T/2π) ~ 5-6, h ≈ 1.1 while mean spacing is ~1 — so grid ≈ zero spacing. d = LT/2π ≈ T·(L/2π) ≈ N. OK.\n3. Build G̃ exactly as E-003 but with GUE points instead of zeta zeros: flat window ψ₀, ϕ = mollified 1_{[−L/2,L/2]}, v_γ = (ϕ̂(γ−α_k)). With L ~ log scale... but for the sine process the \"L\" parameter is free: the limit L→∞ with... The paper's asymptotics: T→∞, L = log(T/2π) → ∞ slowly. The \"limiting spectral distribution\" is presumably the L→∞ limit of the Gabor-Gram eigenvalue distribution over the sine process at critical density h = 2π/L. To simulate: pick moderately large L (say L = 6-12), window of W ≈ d points... this needs points in [T, 2T] with T... hmm, actually we can place the \"window\" anywhere: take sine-process points in [0, T_w] with grid α_k = 2πk/L, k = 0..d−1, d = ⌊LT_w/2π⌋... wait but T here anchors α_0 = T. The construction is shift-invariant; set window [0, W], d = ⌊LW/2π⌋, need W large enough that L = ... no — in the paper L = log(T/2π) is tied to T. For the sine process there's no T; L is a free parameter (the window bandwidth). The limiting distribution as L, W → ∞ with d/W → 1... The paper's HL*(k₀): tr G̃ᵏ = d m_k(1)(1+o(1)) — presumably in the joint limit. So simulate: L = 8, 10, 12; W such that d ~ 1000-2000 → W = 2πd/L.\n\nDirect approach to compute the limiting moments without full simulation: the EXPECTED Gram matrix over the sine process: E[G̃] = Toeplitz with entries ∝ sinc — eigenvalue distribution of THAT (deterministic) is the prolate-spheroidal one. But tr G̃ᵏ has fluctuations; E[tr G̃ᵏ] ≠ tr(E G̃)ᵏ — the k=2 moment includes the variance (that's where 4/3 vs 1 comes from!). Indeed E‖G̃‖² = ‖EG̃‖² + Var-terms; the 4/3 = 1 + 1/3 where 1/3 comes from the sine-kernel |K(x,y)|² term (the \"1−sinc²\" form factor). So the moments m_k are genuinely random-matrix quantities. Monte Carlo with GUE is the way.\n\nSimplest faithful Monte Carlo: use the actual paper construction with zeta zeros themselves but... no, we want sine-process. GUE sampling: n=1000 GUE: eigenvalues, unfold (middle 60%: scale so spacing 1), take a window of W consecutive unfolded eigenvalues ~ 800-1500, build Gram with L = something... but L should grow; with W = 1200 points and L = 10: d = LW/2π ≈ 1910 > W?? Hmm d = LW/(2π): L=10 → d ≈ 1.59W — that's not the paper's d ≈ N. In the paper d = ⌊LT/2π⌋ and N(T,2T) ≈ TL/2π·(1+2log2/L−1/L...) — d ≈ N at leading order because N(T,2T) = (T/2π)(L + 2log2 − 1) + ... hmm wait: N(T,2T) = (T/2π)(log(T/2π) + 2log2) − T/2π + O(log T) = (T/2π)(L − 1 + 2log2) ≈ (T/2π)(L − 0.386). And d = LT/2π. So d/N → L/(L−0.386) ≈ 1.077 at L=5.76 → d ≈ 1.08N. OK.\n\nFor the simulation, parameters: unfolded sine-process points in [0, W] (mean spacing 1, so ~W points), L free ≥ 2-ish, grid h = 2π/L, d = ⌊LW/2π⌋. With W = 600, L = 6: d = 573. G̃ = (aL²)⁻¹ VᵀV, V = ϕ̂(γ_j − α_k) matrix (600×573). Then normalized eigenvalue moments: (1/d)tr G̃ᵏ... wait \"tr G̃ᵏ = d·m_k\" so m_k = (1/d)tr G̃ᵏ = mean of eigenvalues^k. Compute for k=1..4, average over ~50-200 GUE samples, extrapolate L.\n\nExpected: m_1 ≈ 1, m_2 ≈ 4/3, m_3 ≈ 2, m_4 ≈ 13/4.\n\nChristoffel Λ₂(0;1): with moments m_0..m_4 (m_0=1), the Christoffel function at 0 for the moment sequence: Λ₂(0) = 1/(Σ_{i,j≤2} (M⁻¹)_{ij}·0^{i+j})... The Christoffel function Λ_m(x) = 1/Σ_{j=0}^{m} p_j(x)² (orthonormal polys w.r.t. the distribution); at x=0: Λ₂(0) = 1/(p₀(0)² + p₁(0)² + p₂(0)²). Computable from the moment matrix M = (m_{i+j})_{i,j≤2}: Λ₂(0) = (M⁻¹)_{00}·... precisely Λ_m(x) = 1/(u*M⁻¹u) with u = (1, x, x²): at x=0: Λ₂(0) = 1/(M⁻¹)_{00}. With m_0=1, m_1=1, m_2=4/3, m_3=2, m_4=13/4: M = [[1,1,4/3],[1,4/3,2],[4/3,2,13/4]]. Compute (M⁻¹)_{00}: det stuff. M⁻¹_{00} = cofactor/det = det([[4/3,2],[2,13/4]])/det(M) = (4/3·13/4 − 4)/det(M) = (13/3 − 4)/det(M) = (1/3)/det(M). det(M) = 1·(4/3·13/4−4) − 1·(13/4−8/3) + 4/3·(2−16/9) = (1/3) − (13/4−8/3) + (4/3)(2/9) = 1/3 − (39/12−32/12) + 8/27 = 1/3 − 7/12 + 8/27. LCD 108: 36/108 − 63/108 + 32/108 = 5/108. So M⁻¹_00 = (1/3)/(5/108) = 36/5. Λ₂(0) = 5/36 ✔✔ — matches the paper! And the certificate: 1 − Λ₂(0) = 31/36... hmm the paper says HL*(4) gives 13/18 = 26/36. Note 1 − Λ₂(0) = 31/36 ≠ 13/18. Hmm — maybe n₊(G̃)/d ≥ 1 − Λ_m(0) with... 1 − 5/36 = 31/36 ≈ 0.861. But paper says HL*(4) → 13/18 ≈ 0.7222 for simple-on-line. The rank bound gives rank P₁... n₊(G̃)/d ≥ 31/36 would be on-line zeros fraction; simple needs another step (the CMS bound applies to n₊; converting to N^s_0 uses the chain with charges). Whatever — 13/18 = (4 − R(4)... let me see: 13/18 = (4−4/3... no. (2 − R₂)/... With third moment the certificate improves via CMS. Don't need to fully rederive — the numeric reproduction of 5/36 from moments 1,4/3,2,13/4 is already verified by hand above! \n\nSo: Λ₂(0) = 5/36 VERIFIED by hand from the moment matrix (computation above: det = 5/108, cofactor = 1/3, ratio inverse = 5/36 ✔).\n\nAlso 13/18: maybe via the weight ω(m) cubic... §7.3 mentions 13/18? No, §7.2(f): HL*(4) → 13/18 for N^s_0/N. The one-sided CMS bound on n₊... Let me not over-derive; record 5/36 verified; 13/18 derivation path noted for later.\n\nWhat distribution has moments 1, 4/3, 2, 13/4? Try X ~ 1 − cos? Density on [0,2]: ρ(x) = (1/π)... arcsine on [0,2]: ρ = 1/(π√(x(2−x))): mean 1 ✔; m_2 = 3/2 ≠ 4/3. No. Semi-circle on [0,2] (radius 1 centered 1): mean 1, m_2 = 1 + 1/4 = 5/4 ≠ 4/3. No. Uniform[0,2]: 1, 4/3, 2, 16/5: first three match! m_4 = 13/4 = 3.25 vs 16/5 = 3.2 — so close but distinct: NOT uniform. Interesting: matches uniform to 3rd order. The distribution is a perturbation of uniform[0,2].\n\nMonte Carlo should reproduce m_4 = 3.25 to distinguish. Simulation size: need relative accuracy ~1% on m_4. Feasible.\n\nPlan E-004: GUE Monte Carlo:\n- Sample GUE(400), eigenvalues, unfold using exact GUE mean spacing: for bulk eigenvalues, unfolded x_i = λ_i·√(2n)/π·... Simpler: use known GUE semicircle: density ρ(λ) = (1/2π)√(4n−λ²)... for GUE with off-diagonal variance 1, eigenvalues in [−2√n, 2√n]; bulk mean spacing at 0: 1/ρ(0) = 2π/(2√n·... ρ(0) = (1/π)√n·... Let me just take middle eigenvalues and rescale by mean local spacing: spacings of consecutive middle eigenvalues ≈ π/√n·... I'll unfold numerically: fit the cumulative. Standard: take eigenvalues λ_(0.3n..0.7n), rescale by dividing by mean spacing.\n- Window [0,W]: take W_unf consecutive unfolded points ≈ 600.\n- L ∈ {4, 6, 8}: build G̃ as in E-003 (flat window ψ₀; note: the m_k(1) presumably for ψ₀ — the paper says \"sine-kernel Gram matrix\" — flat window corresponds to sinc kernel ✔).\n- Eigenvalues of G̃, normalized moments m_k = mean(λᵏ), k=1..4.\n- Average over 100 samples; report mean ± std; compare to 1, 4/3, 2, 13/4.\n\nSubtlety: with L small (4-8), the mollified ϕ ≈ indicator — the \"sine kernel\" is the Fourier transform of the flat window: ϕ̂(ξ) ∝ sinc(Lξ/2)·L — yes sinc with parameter L/2.\n\nAlso d vs N_points: in window [0,W] with W points and grid d = LW/2π: for L=6: d = 0.955W < W; fine.\n\nPerformance: d ~ 573 at L=6, W=600; V is (W+d')×d — small; 100 samples × 3 L values: trivial.\n\nBut wait — critical check: is the \"sine-kernel Gram matrix\" in the paper built with the smoothed ϕ or the sharp indicator? Sharp indicator ⟹ ϕ̂ = sinc exactly (\"sine kernel\" name suggests sharp). The mollification is a technicality; the limiting distribution should use the sharp cutoff: v entries = sinc((γ−α_k)L/2)... with sharp ϕ = 1_{[−L/2,L/2]}: ϕ̂(ξ) = 2sin(Lξ/2)/ξ = L·sinc(Lξ/2). Then a = ‖ϕ‖²/L = 1. G̃ = (1/L²)VᵀV with V = sinc matrix. I'll offer both but sharp is cleaner: call it the sinc-Gram matrix directly: G̃_{kk'} = (1/L²)Σ_γ sinc(L(γ−α_k)/2)sinc(L(γ−α_k')/2)·L² = Σ_γ sinc(L(γ−α_k)/2)sinc(L(γ−α_k')/2). With sinc(x) := sin(x)/x.\n\nGood—E-004 will compute both ways? Keep simple: sharp sinc (matches \"sine-kernel Gram\" name).\n\nWrite E-004 script, run, compare. Also record the hand-verification of Λ₂(0) = 5/36 in D-010 (short derivation note) and identify density (near-uniform perturbation).\n\nLet me write the script.","signature":""},{"type":"text","text":"Proceeding with **task 1.3 warm-up (thread T1)**. First, a quick hand result before any numerics — the Christoffel value can be verified analytically from the moment sequence:\n\nWith m₀..m₄ = 1, 1, 4/3, 2, 13/4, the moment matrix M = (m_{i+j})_{i,j≤2} has det(M) = 5/108 and (M⁻¹)₀₀ = (1/3)/(5/108) = 36/5, so **Λ₂(0) = 1/(M⁻¹)₀₀ = 5/36 ✔ — matches the paper exactly.** One number already reproduced by hand.\n\nAlso notable: 1, 4/3, 2 matches Uniform[0,2] moments 2^k/(k+1) to third order, but m₄ = 13/4 = 3.25 ≠ 16/5 = 3.2 — the limiting density is a *near-uniform perturbation*. Now the Monte Carlo confirmation (E-004): sample GUE spectra (bulk = sine process), build the sinc-Gram matrix exactly as in §7.2(f), measure the eigenvalue moments."},{"type":"tool_use","id":"Write_188","name":"Write","input":{"file_path":"D:\\Python\\Quantum\\research\\numerics\\scripts\\E-004-sine-moments.py","content":"\"\"\"E-004 — Thread T1 / task 1.3 warm-up: reproduce the limiting spectral\nmoments m_k(1) = 1, 4/3, 2, 13/4 (k <= 4) of the sine-kernel Gram matrix\nover the sine process (paper §7.2(f)), and the Christoffel value\nLambda_2(0) = 5/36.\n\nConstruction (flat window; \"sine-kernel Gram matrix\"):\n  points gamma: sine-process sample (GUE bulk eigenvalues, unfolded to\n  mean spacing 1) in a window of W points;\n  grid alpha_k = 2 pi k / L at critical density h = 2pi/L, k = 0..d-1,\n  d = floor(L W / 2 pi);\n  G_tilde = V^T V with V[j,k] = sinc(L(gamma_j - alpha_k)/2),\n  sinc(x) = sin(x)/x   (sharp cutoff; a = ||phi||^2/L = 1).\nMoments: m_k = (1/d) tr G_tilde^k  = mean of eigenvalues^k.\n\nHand check (already verified, D-010): from m_0..m_4 = 1,1,4/3,2,13/4 the\nmoment matrix has det = 5/108 and (M^-1)_00 = 36/5, so Lambda_2(0) = 5/36.\n\nOutputs: research/numerics/data/E-004-sine-moments.csv (+ .json)\n\"\"\"\nimport json\nimport sys\nimport time\nfrom datetime import datetime, timezone\nfrom pathlib import Path\n\nimport numpy as np\n\nREPO = Path(__file__).resolve().parents[3]\nDATA = REPO / \"research\" / \"numerics\" / \"data\"\nCSV = DATA / \"E-004-sine-moments.csv\"\nJSON = DATA / \"E-004-sine-moments.json\"\n\nRNG_SEED = 20260817\nN_GUE = 400          # GUE size per sample (bulk used)\nW = 600              # window length in unfolded points\nL_VALUES = (4.0, 6.0, 8.0)\nN_SAMPLES = 60\nTARGET = {1: 1.0, 2: 4.0 / 3.0, 3: 2.0, 4: 13.0 / 4.0}\n\n\ndef gue_bulk_points(rng) -> np.ndarray:\n    \"\"\"Unfolded GUE bulk eigenvalues (mean spacing 1), W consecutive.\"\"\"\n    A = rng.standard_normal((N_GUE, N_GUE)) + 1j * rng.standard_normal((N_GUE, N_GUE))\n    H = (A + A.conj().T) / 2.0\n    ev = np.linalg.eigvalsh(H)  # support ~ [-2 sqrt(2n), 2 sqrt(2n)] for this normalization? no:\n    # entries variance 1 (complex): semicircle radius 2*sqrt(n)\n    ev = ev / np.sqrt(N_GUE)    # now support [-2, 2], bulk density (1/pi)... unfold:\n    # semicircle density rho(x) = (1/(2 pi)) sqrt(4 - x^2); unfold via CDF\n    # take middle 60% and rescale by local mean spacing at 0: spacing(0) = 1/rho(0) = pi\n    lo, hi = int(0.2 * N_GUE), int(0.8 * N_GUE)\n    bulk = ev[lo:hi]\n    bulk = bulk * (1.0 / np.pi) * N_GUE / (0.6 * N_GUE) * len(bulk) * (np.pi / len(bulk))\n    # simpler robust unfolding: rescale so mean consecutive spacing is exactly 1\n    bulk = bulk / np.mean(np.diff(bulk))\n    mid = len(bulk) // 2\n    return bulk[mid - W // 2: mid + W // 2]\n\n\ndef sinc(x):\n    return np.sinc(x / np.pi)  # np.sinc(u) = sin(pi u)/(pi u) -> sin(x)/x\n\n\ndef main() -> int:\n    if CSV.exists() and \"--force\" not in sys.argv:\n        print(f\"cached: {CSV} exists (use --force to recompute)\")\n        return 0\n    rng = np.random.default_rng(RNG_SEED)\n    t0 = time.time()\n\n    # accumulate moments per L\n    acc = {L: {k: [] for k in (1, 2, 3, 4)} for L in L_VALUES}\n    for s in range(N_SAMPLES):\n        pts = gue_bulk_points(rng)\n        pts = pts - pts[0]\n        Wpts = pts[-1]\n        for L in L_VALUES:\n            h = 2 * np.pi / L\n            d = int(L * Wpts / (2 * np.pi))\n            alpha = h * np.arange(d)\n            Z = pts[:, None] - alpha[None, :]\n            V = sinc(L * Z / 2.0)\n            G = (V.T @ V) / L**2\n            ev = np.linalg.eigvalsh(G)\n            for k in (1, 2, 3, 4):\n                acc[L][k].append(float(np.mean(ev**k)))\n        if (s + 1) % 20 == 0:\n            print(f\"  {s+1}/{N_SAMPLES} samples ({time.time()-t0:.0f}s)\",\n                  flush=True)\n\n    rows = []\n    checks = {}\n    for L in L_VALUES:\n        row = {\"L\": L}\n        for k in (1, 2, 3, 4):\n            arr = np.array(acc[L][k])\n            m, sd = arr.mean(), arr.std()\n            row[f\"m{k}\"] = f\"{m:.5f}\"\n            row[f\"m{k}_std\"] = f\"{sd:.5f}\"\n            row[f\"m{k}_target\"] = f\"{TARGET[k]:.5f}\"\n            row[f\"m{k}_relerr\"] = f\"{abs(m / TARGET[k] - 1):.5f}\"\n        rows.append(row)\n        # generous tolerances: finite-L + finite-W + sampling effects\n        checks[f\"L={L}\"] = {\n            f\"m{k}_within_6pct\": bool(abs(np.mean(acc[L][k]) / TARGET[k] - 1) < 0.06)\n            for k in (1, 2, 3, 4)\n        }\n        print(f\"L={L}: \" + \"  \".join(\n            f\"m{k}={np.mean(acc[L][k]):.4f}+-{np.std(acc[L][k]):.4f}\"\n            f\" (target {TARGET[k]:.4f})\" for k in (1, 2, 3, 4)), flush=True)\n\n    # Christoffel from simulated moments (L=8): Lambda_2(0) = 1/(M^-1)_00\n    Lmax = L_VALUES[-1]\n    m = [1.0] + [float(np.mean(acc[Lmax][k])) for k in (1, 2, 3, 4)]\n    M = np.array([[m[i + j] for j in range(3)] for i in range(3)])\n    lam2 = 1.0 / np.linalg.inv(M)[0, 0]\n    print(f\"Christoffel Lambda_2(0) from simulated moments (L={Lmax}): \"\n          f\"{lam2:.5f}  (target 5/36 = {5/36:.5f})\")\n    checks[\"christoffel\"] = {\"within_10pct\": bool(abs(lam2 / (5 / 36) - 1) < 0.10)}\n\n    DATA.mkdir(parents=True, exist_ok=True)\n    with CSV.open(\"w\", encoding=\"utf-8\", newline=\"\") as f:\n        f.write(\",\".join(rows[0].keys()) + \"\\n\")\n        for r in rows:\n            f.write(\",\".join(str(v) for v in r.values()) + \"\\n\")\n    meta = {\n        \"experiment\": \"E-004\",\n        \"task\": \"1.3 warm-up (thread T1)\",\n        \"description\": \"Monte Carlo reproduction of sine-kernel Gram matrix \"\n                       \"spectral moments m_k(1)=1,4/3,2,13/4 and Christoffel \"\n                       \"Lambda_2(0)=5/36 (paper §7.2(f))\",\n        \"script\": \"research/numerics/scripts/E-004-sine-moments.py\",\n        \"parameters\": {\"N_GUE\": N_GUE, \"W\": W, \"L_values\": list(L_VALUES),\n                       \"n_samples\": N_SAMPLES, \"seed\": RNG_SEED},\n        \"christoffel_simulated\": lam2,\n        \"checks\": checks,\n        \"elapsed_seconds\": round(time.time() - t0, 1),\n        \"timestamp_utc\": datetime.now(timezone.utc).isoformat(),\n        \"row_count\": len(rows),\n    }\n    JSON.write_text(json.dumps(meta, indent=2), encoding=\"utf-8\")\n    ok = all(all(c.values()) for c in checks.values())\n    print(f\"wrote {CSV.name}, {JSON.name} — \"\n          f\"{'ALL CHECKS PASS' if ok else 'CHECK FAILURE (inspect)'}\")\n    return 0 if ok else 1\n\n\nif __name__ == \"__main__\":\n    raise SystemExit(main())\n"},"extra_content":null}]