← Session 5 — Phase 0: paper deep-read + Bombieri–Lagarias literature
thinking (7718 chars)
Now write E-009a script. Plan:
```python
"""E-009a — ..."""
```
Key functions:
- von Mangoldt sieve up to N = ceil(e^10) ≈ 22027 → use 22100.
- prime_sum(t) via cumsum arrays (float64) and mpmath exact version for small t (direct loop).
- g_direct(t): mpmath high precision:
-4*(exp(t/2)+exp(-t/2)-2) + primesum(|t|) - (|t|/2)*(digamma(0.25)-log(pi)) - 0.25*(lerchphi(1,2,0.25) - exp(-|t|/2)*lerchphi(exp(-2*|t|),2,0.25))
Note prime sum Σ_{n≤e^{|t|}} Λ(n)/√n (|t|−log n).
mpmath lerchphi signature: lerchphi(z, s, a). ✔
Validations (dps 30):
1. lerch convention: lerchphi(1,2,1/4) == polygamma(1, 1/4).
2. Evenness: g(t)−g(−t) = 0 for t in {0.7, 2.3, 5.1}.
3. A-limit: g(t)/t − ½log t → A for t = 1e-2, 1e-3, 1e-4 (t below ½log2 so prime sum = 0). Error should be O(t): check |err(1e-4)| < 1e-3 and decreasing ~10x each step.
4. r-ratio: r(t)/t² convergence for t = 0.25, 0.125, ..., ~2e-3: differences decrease; report.
Tabulation float64: t in [0,10], h=1e-3, columns: t, g(t), prime_part(t), smooth check maybe. Columns: t,g,prime_sum. Compute g via float64 path: use math/cumsum; lerchphi float64 via mpmath at dps 15 — 10k lerchphi evals might be slowish but OK (~ each lerchphi with small z converges fast; maybe 1ms each → 10s). Fine. Actually we can vectorize the Lerch sum directly: Φ(e^{−2t},2,1/4) = Σ_{k=0}^{K} e^{−2tk}/(k+.25)², converges like e^{−2tK}; for tiny t need many terms (for t=1e-3, need K ~ 1e4 for tail e^{−2tK}<1e-12 → K ~ 14000). Hmm, for small t the direct sum is slow-converging. Use mpmath lerchphi which handles it. 10001 evals — test speed; maybe a few minutes. Acceptable. Alternatively precompute on a coarse mesh + interpolate... no, keep exact.
Actually, better: compute the float64 table in numpy vectorized lerch for t ≥ some cutoff, and mpmath for t < cutoff. Meh — just use mpmath for all; if too slow, adjust. Let me estimate: lerchphi(z,s,a) for z close to 1 uses... mpmath may be slow (it does series acceleration, maybe 0.1-1ms). 10k × 1ms = 10s. OK.
For float64 speed, also possible: note the combination −¼[Φ(1,2,¼) − e^{−t/2}Φ(e^{−2t},2,¼)] is small for small t (cancellation!). At t=1e-3, Φ(e^{−2t},2,¼) ≈ Φ(1,2,¼) − 2tΦ(1,1,¼)... Φ(1,1,¼) diverges, so difference ~ t log t·(−?) Anyway cancellation of leading digits: Φ values ~17, difference ~ 1e-3·log ~ 1e-3 → lose ~4 digits at t=1e-3, worse at smaller t. float64 gives ~12 reliable digits at t=1e-3, at t=1e-4 ~11... For Galerkin purposes that's plenty. But for CSV table, note the caveat. Alternatively store the combination computed at dps 30 then converted to float — better. Use mpmath dps 30 for the whole table computation? 10k evals with prime cumsum fast — the lerchphi at dps 30 maybe 10ms each → 100s. Acceptable but let's do dps 30 for t < 0.1 (100 points) and dps 20 elsewhere. Actually simpler: do everything at dps 25; prime sum via float cumsum (error negligible relative). Let's estimate lerchphi cost... I'll just write it with a progress print and run; if slow, optimize.
Hmm — actually a smarter approach for the table: the singular structure is at t=0 only. We can precompute the Lerch-combination on the mesh efficiently by noting Φ(e^{−2t},2,1/4) at mesh points t_k = k h: it's Σ_k z_k^m/(m+1/4)² with z = e^{−2h}. That's a power series evaluated at geometric points — could use fast multipoint, overkill. Just mpmath.
Let me now write the script.
Also A constant: A = 0.5*(log(2*pi) + euler - 1) = 0.70754636568...? compute: 1.83787706641 + 0.57721566490 = 2.41509273131; −1 = 1.41509273131; /2 = 0.70754636565. D-018 says ≈ 0.707546 ✔.
digamma(0.25): −γ − π/2·cot(π/4) − 3log2 = −0.5772156649 − 1.5707963268 − 2.07944154168 = −4.22745353338. −log π shift: ψ(1/4) − logπ = −4.22745353338 − 1.14472988585 = −5.37218341923. Coefficient −(|t|/2)(...) = +2.68609170962·|t|.
Write CSV with columns: t, g_t, prime_sum_t. Plus maybe g minus singular parts? Keep simple: t,g,prime_sum.
Checks:
- check_lerch_identity
- check_evenness (max diff < 1e-25)
- check_A_limit (err at 1e-4 < 2e-4 and monotone decreasing ~ ×10)
- check_r_ratio_convergence (|last two ratios| relative diff < 0.02)
- check_table_finite & monotone-ish? Table sanity: g(10) finite, no NaN.
JSON sidecar with key values: A constant, r2 = lim r(t)/t² estimate, sample g values, digamma coefficient.
Threshold for A-limit test: err(t) ≈ C·t with C = −½r''(0)·... err = (A + (r(t)/t)) − A... g(t)/t − ½ log t = A + r(t)/t = A + O(t). At t=1e-2: O(1e-2·c). We don't know c = ½r''(0) magnitude; from the −t² contribution alone c≈−1 → err ~ 1e-2 at t=1e-2? Wait r(t)/t = O(t), so at 1e-2 err ~ c·1e-2 with c~1 → 1e-2. Hmm, that's a big tolerance. Check ratio err(1e-3)/err(1e-2) ≈ 0.1. Test: |err(1e-4)| < 3e-3 and err ratios in [0.05, 0.2]. Empirically the r-ratio test will pin c: lim r(t)/t² = c, then expected err(1e-4) = c·1e-4. I'll compute c first, then test A-limit err ≈ c·t within 20%.
Order: compute r-ratio first to get c, then A-limit test against A + c·t.
Actually r(t)/t² → c, and g(t)/t − ½log t − A = r(t)/t = c·t + o(t). So test |observed − c·t| / |c·t| small. Good.
But note: the A-limit test uses t below ½ log 2 ≈ 0.3466; 1e-2..1e-4 fine. r-ratio uses t up to 0.25 (also below 0.3466 ✔).
Now the table: use mpmath mp.dps = 25 for lerchphi, but prime sum from float cumsum (error ~1e-13·magnitude, fine at 25dps? It contaminates the 25-digit value — but we only print ~17 digits to CSV. Set dps 20 for speed.)
Runtime estimate for 10001 lerchphi calls at dps 20: unknown; add time printout. If too slow I'll subsample near 0. Let me include a fallback: compute lerch combination via direct series when e^{−2t} < 0.5 (t > 0.35): fast numpy vectorized. For t ≤ 0.35 (350 points) use mpmath. Actually z = e^{−2t} < 0.5 means terms decay by ≥0.5^k — need only ~50 terms for 1e-15. For t>0.35 vectorized numpy over 10k points × 60 terms trivial. For t ≤ 0.35: 351 mpmath calls — fast.
Also at t=0: g=0 exactly, set directly.
Numpy direct series: Phi(z) = Σ_{k=0}^{K} z^k/(k+0.25)^2. For array of z values (all < 0.5), K=80.
Then combination C(t) = −0.25*(Phi1 − exp(−t/2)*Phi(exp(−2t))) where Phi1 = polygamma(1,0.25) float64 = 17.19732915450756... get via mpmath once at dps 30, hardcode as computed value in JSON.
Also −4(cosh term): −4*(exp(t/2)+exp(−t/2)−2).
digamma term float: coeff = 2.6860917096128... compute via mpmath once.
Prime sum via cumsum: logs = log(n) for n=2..N; sieve Λ; lam = Λ/√n; C1 = cumsum(lam), C0 = cumsum(lam*log n); m = np.searchsorted(logs, t, 'right'); P = C1[m−1]*t − C0[m−1] (if m>0).
For high-precision validation prime sum: direct mpmath loop over prime powers ≤ e^t (t ≤ 0.25 → n ≤ 1, wait e^0.25 = 1.28 < 2, prime sum = 0!). Oh nice — for ALL validation t ≤ 0.3466 the prime sum vanishes. So validation needs no prime sum at all. But evenness test at t=2.3, 5.1 does need it — use mpmath loop there (n ≤ e^5.1 = 164, small).
Sieve: N = int(exp(10))+2 = 22028. lam array.
Now write the script.
def von_mangoldt_sieve(N): returns array lam_log[0..N] where lam_log[n] = log p if n = p^k.
```python
lam_log = np.zeros(N+1)
is_composite = np.zeros(N+1, bool)
for p in range(2, N+1):
if not is_composite[p]:
pk = p
while pk <= N:
lam_log[pk] = np.log(p)
# mark multiples of p composite? standard sieve: mark p-multiples
pk *= p
for m in range(2*p, N+1, p):
is_composite[m] = True
```
Wait marking must happen for all p; do standard sieve first for primes, then powers.
CSV rows 10001 — fine.
JSON checks and key values.
Let me also add g evaluation at interesting points for the record: g(½log2) (first prime kink), g(log 2), g(log 4)... record in JSON.
Write the file now.
Conventions clear. Writing the E-009a script now.
tool: Write
input
{
"file_path": "D:\\Python\\Quantum\\research\\numerics\\scripts\\E-009-suzuki-screw-function.py",
"content": "\"\"\"E-009a \u2014 Task 3.4c (Q-3.4c, D-018 plan step 1): Suzuki's screw function g(t).\n\nSuzuki (arXiv:2606.09096v1, eq. 1.3) \u2014 explicit screw function of zeta:\n\n g(t) = -4(e^{t/2} + e^{-t/2} - 2)\n + Sum_{n <= e^{|t|}} (Lambda(n)/sqrt(n)) (|t| - log n)\n - (|t|/2)(psi(1/4) - log pi)\n - (1/4)[ Phi(1,2,1/4) - e^{-|t|/2} Phi(e^{-2|t|}, 2, 1/4) ]\n\n(psi = digamma, Phi = Hurwitz-Lerch transcendent Phi(z,s,a) =\nSum_{k>=0} z^k/(k+a)^s; Lambda = von Mangoldt; prime sum = 0 for\n|t| < (1/2) log 2.)\n\nNear-origin expansion (Suzuki eq. 2.2):\n g(t) = (1/2)|t| log|t| + A|t| + (prime sum) + r(t),\n r even C^2, r = O(t^2),\n A = (1/2)(log(2pi) - psi(2)) = (1/2)(log(2pi) + gamma_E - 1)\n = 0.7075463656...\n\ng is a Krein-Langer screw function IFF RH (this is the Cor 1.6 bridge).\n\nPart A (validation, mpmath dps 30):\n A1. Lerch convention: Phi(1,2,1/4) == psi'(1/4) = polygamma(1, 1/4).\n A2. Evenness: g(-t) == g(t) at t = 0.7, 2.3, 5.1 (exercises prime sum).\n A3. A-limit: g(t)/t - (1/2) log t -> A as t -> 0+ (prime sum = 0 there);\n residual must equal r(t)/t = c2*t + o(t) with c2 from A4.\n A4. r-curvature: r(t)/t^2 -> c2 (constant) as t -> 0+, t <= 0.25.\n\nPart B (tabulation, float64 + mpmath hybrid):\n g(t) on t in [0, 10], step 1e-3 (mesh for the E-009b Galerkin stage;\n g is even so t >= 0 suffices). Prime-sum part via exact cumsum arrays\n (P(t) piecewise linear in t with kinks at log n). Lerch combination:\n direct numpy series for t > 0.35 (z = e^{-2t} < 0.5, fast), mpmath\n lerchphi for t <= 0.35 (slow-converging near z = 1).\n\nOutputs: research/numerics/data/E-009-suzuki-screw-function.csv (+ .json)\n\"\"\"\nimport json\nimport sys\nimport time\nfrom datetime import datetime, timezone\nfrom pathlib import Path\n\nimport numpy as np\nimport mpmath as mp\n\nREPO = Path(__file__).resolve().parents[3]\nDATA = REPO / \"research\" / \"numerics\" / \"data\"\nCSV = DATA / \"E-009-suzuki-screw-function.csv\"\nJSON = DATA / \"E-009-suzuki-screw-function.json\"\n\nT_MAX = 10.0 # tabulation range [0, T_MAX] (g even)\nH = 1e-3 # mesh step\nDPS = 30 # validation precision\nA_CONST = 0.5 * (np.log(2 * np.pi) + np.euler_gamma - 1.0) # eq. 2.2 constant\n\n\n# ---------------------------------------------------------------- von Mangoldt\ndef mangoldt_sieve(n_max: int) -> np.ndarray:\n \"\"\"lam_log[n] = log p if n = p^k (prime power), else 0.\"\"\"\n lam_log = np.zeros(n_max + 1)\n is_comp = np.zeros(n_max + 1, dtype=bool)\n for p in range(2, n_max + 1):\n if not is_comp[p]:\n if p * p <= n_max:\n is_comp[p * p :: p] = True\n pk = p\n while pk <= n_max:\n lam_log[pk] = np.log(p)\n pk *= p\n return lam_log\n\n\nN_MAX = int(np.exp(T_MAX)) + 2\nLAM_LOG = mangoldt_sieve(N_MAX)\nNS = np.arange(2, N_MAX + 1)\nLOGS = np.log(NS)\nLAM = LAM_LOG[2:] / np.sqrt(NS) # Lambda(n)/sqrt(n)\nC1 = np.cumsum(LAM) # Sum_{n<=x} Lambda(n)/sqrt(n)\nC0 = np.cumsum(LAM * LOGS) # Sum_{n<=x} Lambda(n) log(n)/sqrt(n)\n\n\ndef prime_sum_vec(t: np.ndarray) -> np.ndarray:\n \"\"\"P(t) = Sum_{n <= e^t} (Lambda(n)/sqrt(n))(t - log n), vectorized.\"\"\"\n m = np.searchsorted(LOGS, t, side=\"right\") # count of n with log n <= t\n out = np.zeros_like(t)\n pos = m > 0\n idx = m[pos] - 1\n out[pos] = C1[idx] * t[pos] - C0[idx]\n return out\n\n\ndef prime_sum_mp(t: mp.mpf) -> mp.mpf:\n s = mp.mpf(0)\n n = 2\n while mp.log(n) <= t:\n s += (LAM_LOG[n] / mp.sqrt(n)) * (t - mp.log(n))\n n += 1\n return s\n\n\n# ------------------------------------------------------------------- g(t), mp\ndef g_mp(t) -> mp.mpf:\n t = mp.mpf(t)\n a = abs(t)\n s = -4 * (mp.exp(a / 2) + mp.exp(-a / 2) - 2)\n s += prime_sum_mp(a)\n s -= (a / 2) * (mp.digamma(mp.mpf(1) / 4) - mp.log(mp.pi))\n s -= mp.mpf(1) / 4 * (mp.lerchphi(1, 2, mp.mpf(1) / 4)\n - mp.exp(-a / 2) * mp.lerchphi(mp.exp(-2 * a), 2,\n mp.mpf(1) / 4))\n return s\n\n\n# --------------------------------------------------------------- g(t), float\nKK = np.arange(0, 80)\nDENOM = (KK + 0.25) ** 2\n\n\ndef lerch_phi2_fast(z: np.ndarray) -> np.ndarray:\n \"\"\"Phi(z,2,1/4) by direct series; valid for |z| <= 0.5 (t >= 0.35).\"\"\"\n powers = np.power.outer(z, KK) # (n_z, K)\n return (powers / DENOM).sum(axis=1)\n\n\ndef g_float_table(t: np.ndarray) -> np.ndarray:\n g = -4.0 * (np.exp(t / 2) + np.exp(-t / 2) - 2.0)\n g += prime_sum_vec(t)\n dig_coeff = -0.5 * float(mp.digamma(mp.mpf(1) / 4) - mp.log(mp.pi))\n g += dig_coeff * t\n phi1 = float(mp.lerchphi(1, 2, mp.mpf(1) / 4)) # = polygamma(1, 1/4)\n lerch_comb = np.empty_like(t)\n slow = t <= 0.35\n if slow.any():\n mp.dps = 25\n for i in np.flatnonzero(slow):\n ti = mp.mpf(float(t[i]))\n val = (mp.lerchphi(1, 2, mp.mpf(1) / 4)\n - mp.exp(-ti / 2) * mp.lerchphi(mp.exp(-2 * ti), 2,\n mp.mpf(1) / 4))\n lerch_comb[i] = float(val) / -4.0\n mp.dps = DPS\n if (~slow).any():\n z = np.exp(-2.0 * t[~slow])\n lerch_comb[~slow] = -0.25 * (phi1 - np.exp(-t[~slow] / 2)\n * lerch_phi2_fast(z))\n g += lerch_comb\n g[t == 0.0] = 0.0\n return g\n\n\n# ------------------------------------------------------------------- main\ndef main() -> int:\n if CSV.exists() and \"--force\" not in sys.argv:\n print(f\"cached: {CSV} exists (use --force to recompute)\"); return 0\n t0 = time.time()\n mp.dps = DPS\n\n # ---- A1: Lerch convention\n phi1_mp = mp.lerchphi(1, 2, mp.mpf(1) / 4)\n trig1_mp = mp.polygamma(1, mp.mpf(1) / 4)\n a1_err = abs(phi1_mp - trig1_mp)\n check_a1 = bool(a1_err < mp.mpf(\"1e-25\"))\n\n # ---- A2: evenness (t values exercising the prime sum)\n ev_ts = [mp.mpf(\"0.7\"), mp.mpf(\"2.3\"), mp.mpf(\"5.1\")]\n ev_err = max(abs(g_mp(t) - g_mp(-t)) for t in ev_ts)\n check_a2 = bool(ev_err < mp.mpf(\"1e-20\"))\n\n # ---- A4: r(t)/t^2 -> c2 (t <= 0.25: prime sum = 0)\n A_mp = (mp.log(2 * mp.pi) + mp.euler - 1) / 2\n r_ts = [mp.mpf(\"0.25\"), mp.mpf(\"0.125\"), mp.mpf(\"0.0625\"),\n mp.mpf(\"0.03125\"), mp.mpf(\"0.015625\")]\n r_ratios = []\n for t in r_ts:\n r = g_mp(t) - (t / 2) * mp.log(t) - A_mp * t # prime sum = 0 here\n r_ratios.append(r / t**2)\n c2 = r_ratios[-1]\n r_gap_last = abs(r_ratios[-1] - r_ratios[-2])\n r_gap_prev = abs(r_ratios[-2] - r_ratios[-3])\n check_a4 = bool(r_gap_last < r_gap_prev) # gaps shrinking => converging\n\n # ---- A3: g(t)/t - (1/2) log t -> A; residual ~ c2 * t\n a3_ts = [mp.mpf(\"1e-2\"), mp.mpf(\"1e-3\"), mp.mpf(\"1e-4\")]\n a3_resid = [g_mp(t) / t - mp.log(t) / 2 - A_mp for t in a3_ts]\n # expected residual ~ c2*t; test relative agreement within 25%\n a3_rel = [abs((res - c2 * t) / (c2 * t)) if c2 != 0 else abs(res)\n for res, t in zip(a3_resid, a3_ts)]\n check_a3 = bool(a3_rel[-1] < 0.25 and abs(a3_resid[-1]) < mp.mpf(\"3e-3\"))\n\n print(f\"A1 lerch==polygamma(1,1/4): err {mp.nstr(a1_err, 5)}\")\n print(f\"A2 evenness max err: {mp.nstr(ev_err, 5)}\")\n print(f\"A4 r(t)/t^2 -> c2 = {mp.nstr(c2, 12)} \"\n f\"(gaps {mp.nstr(r_gap_prev, 3)} -> {mp.nstr(r_gap_last, 3)})\")\n print(f\"A3 residuals g(t)/t - (1/2)log t - A at 1e-2/1e-3/1e-4: \"\n f\"{[mp.nstr(x, 6) for x in a3_resid]} (c2*t predicted)\")\n\n # ---- Part B: table on [0, T_MAX]\n t_mesh = np.arange(0.0, T_MAX + 0.5 * H, H)\n g_vals = g_float_table(t_mesh)\n p_vals = prime_sum_vec(t_mesh)\n check_b = bool(np.all(np.isfinite(g_vals)))\n # spot check: float table vs mpmath at a few points\n spot_ts = [0.123, 1.0, 3.7, 7.0, 10.0]\n spot_err = max(abs(g_mp(x) - mp.mpf(float(g_vals[int(round(x / H))])))\n for x in spot_ts)\n check_spot = bool(spot_err < mp.mpf(\"1e-9\"))\n print(f\"Part B: {len(t_mesh)} points, spot-check err vs mpmath: \"\n f\"{mp.nstr(spot_err, 3)}\")\n\n checks = {\n \"A1_lerch_convention\": check_a1,\n \"A2_evenness\": check_a2,\n \"A3_A_limit\": check_a3,\n \"A4_r_curvature_converges\": check_a4,\n \"B_table_finite\": check_b,\n \"B_spot_check_vs_mpmath\": check_spot,\n }\n\n DATA.mkdir(parents=True, exist_ok=True)\n with CSV.open(\"w\", encoding=\"utf-8\", newline=\"\") as f:\n f.write(\"t,g_t,prime_sum_t\\n\")\n for ti, gi, pi in zip(t_mesh, g_vals, p_vals):\n f.write(f\"{ti:.4f},{gi:.12e},{pi:.12e}\\n\")\n\n meta = {\n \"experiment\": \"E-009a\",\n \"task\": \"3.4c (D-018 step 1): Suzuki screw function g(t)\",\n \"description\": \"mpmath implementation of Suzuki's explicit screw \"\n \"function g(t) (arXiv:2606.09096v1 eq. 1.3); \"\n \"validated against eq. 2.2 near-origin expansion; \"\n \"tabulated on [0,10] step 1e-3 for E-009b Galerkin\",\n \"script\": \"research/numerics/scripts/E-009-suzuki-screw-function.py\",\n \"parameters\": {\"t_max\": T_MAX, \"step\": H, \"validation_dps\": DPS,\n \"n_sieve_max\": N_MAX},\n \"checks\": checks,\n \"key_values\": {\n \"A_eq2.2\": float(A_CONST),\n \"c2_lim_r_over_t2\": float(c2),\n \"digamma_coeff_-half(psi(1/4)-log pi)\":\n float(-(mp.digamma(mp.mpf(1) / 4) - mp.log(mp.pi)) / 2),\n \"Phi(1,2,1/4)\": float(phi1_mp),\n \"g(log2)\": float(g_mp(mp.log(2))),\n \"g(1)\": float(g_mp(1)),\n \"g(5)\": float(g_mp(5)),\n \"g(10)\": float(g_vals[-1]),\n },\n \"caveats\": [\n \"float64 table: Lerch combination has cancellation for small t \"\n \"(~4 digits lost at t=1e-3); spot-checked vs mpmath to 1e-9\",\n \"prime sum is piecewise linear in t (kinks at log n) \u2014 exact \"\n \"via cumsum, no interpolation error\",\n \"g even: only t >= 0 tabulated\",\n ],\n \"elapsed_seconds\": round(time.time() - t0, 1),\n \"timestamp_utc\": datetime.now(timezone.utc).isoformat(),\n \"row_count\": int(len(t_mesh)),\n }\n JSON.write_text(json.dumps(meta, indent=2), encoding=\"utf-8\")\n ok = all(checks.values())\n print(f\"wrote {CSV.name}, {JSON.name} \u2014 \"\n f\"{'ALL CHECKS PASS' if ok else 'CHECK FAILURE: ' + str(checks)}\")\n return 0 if ok else 1\n\n\nif __name__ == \"__main__\":\n raise SystemExit(main())\n"
}
raw content_json (database row)
[{"type":"thinking","thinking":"Now write E-009a script. Plan:\n\n```python\n\"\"\"E-009a — ...\"\"\"\n```\n\nKey functions:\n- von Mangoldt sieve up to N = ceil(e^10) ≈ 22027 → use 22100.\n- prime_sum(t) via cumsum arrays (float64) and mpmath exact version for small t (direct loop).\n- g_direct(t): mpmath high precision: \n -4*(exp(t/2)+exp(-t/2)-2) + primesum(|t|) - (|t|/2)*(digamma(0.25)-log(pi)) - 0.25*(lerchphi(1,2,0.25) - exp(-|t|/2)*lerchphi(exp(-2*|t|),2,0.25))\n\nNote prime sum Σ_{n≤e^{|t|}} Λ(n)/√n (|t|−log n).\n\nmpmath lerchphi signature: lerchphi(z, s, a). ✔\n\nValidations (dps 30):\n1. lerch convention: lerchphi(1,2,1/4) == polygamma(1, 1/4).\n2. Evenness: g(t)−g(−t) = 0 for t in {0.7, 2.3, 5.1}.\n3. A-limit: g(t)/t − ½log t → A for t = 1e-2, 1e-3, 1e-4 (t below ½log2 so prime sum = 0). Error should be O(t): check |err(1e-4)| < 1e-3 and decreasing ~10x each step.\n4. r-ratio: r(t)/t² convergence for t = 0.25, 0.125, ..., ~2e-3: differences decrease; report.\n\nTabulation float64: t in [0,10], h=1e-3, columns: t, g(t), prime_part(t), smooth check maybe. Columns: t,g,prime_sum. Compute g via float64 path: use math/cumsum; lerchphi float64 via mpmath at dps 15 — 10k lerchphi evals might be slowish but OK (~ each lerchphi with small z converges fast; maybe 1ms each → 10s). Fine. Actually we can vectorize the Lerch sum directly: Φ(e^{−2t},2,1/4) = Σ_{k=0}^{K} e^{−2tk}/(k+.25)², converges like e^{−2tK}; for tiny t need many terms (for t=1e-3, need K ~ 1e4 for tail e^{−2tK}<1e-12 → K ~ 14000). Hmm, for small t the direct sum is slow-converging. Use mpmath lerchphi which handles it. 10001 evals — test speed; maybe a few minutes. Acceptable. Alternatively precompute on a coarse mesh + interpolate... no, keep exact.\n\nActually, better: compute the float64 table in numpy vectorized lerch for t ≥ some cutoff, and mpmath for t < cutoff. Meh — just use mpmath for all; if too slow, adjust. Let me estimate: lerchphi(z,s,a) for z close to 1 uses... mpmath may be slow (it does series acceleration, maybe 0.1-1ms). 10k × 1ms = 10s. OK.\n\nFor float64 speed, also possible: note the combination −¼[Φ(1,2,¼) − e^{−t/2}Φ(e^{−2t},2,¼)] is small for small t (cancellation!). At t=1e-3, Φ(e^{−2t},2,¼) ≈ Φ(1,2,¼) − 2tΦ(1,1,¼)... Φ(1,1,¼) diverges, so difference ~ t log t·(−?) Anyway cancellation of leading digits: Φ values ~17, difference ~ 1e-3·log ~ 1e-3 → lose ~4 digits at t=1e-3, worse at smaller t. float64 gives ~12 reliable digits at t=1e-3, at t=1e-4 ~11... For Galerkin purposes that's plenty. But for CSV table, note the caveat. Alternatively store the combination computed at dps 30 then converted to float — better. Use mpmath dps 30 for the whole table computation? 10k evals with prime cumsum fast — the lerchphi at dps 30 maybe 10ms each → 100s. Acceptable but let's do dps 30 for t < 0.1 (100 points) and dps 20 elsewhere. Actually simpler: do everything at dps 25; prime sum via float cumsum (error negligible relative). Let's estimate lerchphi cost... I'll just write it with a progress print and run; if slow, optimize.\n\nHmm — actually a smarter approach for the table: the singular structure is at t=0 only. We can precompute the Lerch-combination on the mesh efficiently by noting Φ(e^{−2t},2,1/4) at mesh points t_k = k h: it's Σ_k z_k^m/(m+1/4)² with z = e^{−2h}. That's a power series evaluated at geometric points — could use fast multipoint, overkill. Just mpmath.\n\nLet me now write the script.\n\nAlso A constant: A = 0.5*(log(2*pi) + euler - 1) = 0.70754636568...? compute: 1.83787706641 + 0.57721566490 = 2.41509273131; −1 = 1.41509273131; /2 = 0.70754636565. D-018 says ≈ 0.707546 ✔.\n\ndigamma(0.25): −γ − π/2·cot(π/4) − 3log2 = −0.5772156649 − 1.5707963268 − 2.07944154168 = −4.22745353338. −log π shift: ψ(1/4) − logπ = −4.22745353338 − 1.14472988585 = −5.37218341923. Coefficient −(|t|/2)(...) = +2.68609170962·|t|.\n\nWrite CSV with columns: t, g_t, prime_sum_t. Plus maybe g minus singular parts? Keep simple: t,g,prime_sum.\n\nChecks:\n- check_lerch_identity\n- check_evenness (max diff < 1e-25)\n- check_A_limit (err at 1e-4 < 2e-4 and monotone decreasing ~ ×10)\n- check_r_ratio_convergence (|last two ratios| relative diff < 0.02)\n- check_table_finite & monotone-ish? Table sanity: g(10) finite, no NaN.\n\nJSON sidecar with key values: A constant, r2 = lim r(t)/t² estimate, sample g values, digamma coefficient.\n\nThreshold for A-limit test: err(t) ≈ C·t with C = −½r''(0)·... err = (A + (r(t)/t)) − A... g(t)/t − ½ log t = A + r(t)/t = A + O(t). At t=1e-2: O(1e-2·c). We don't know c = ½r''(0) magnitude; from the −t² contribution alone c≈−1 → err ~ 1e-2 at t=1e-2? Wait r(t)/t = O(t), so at 1e-2 err ~ c·1e-2 with c~1 → 1e-2. Hmm, that's a big tolerance. Check ratio err(1e-3)/err(1e-2) ≈ 0.1. Test: |err(1e-4)| < 3e-3 and err ratios in [0.05, 0.2]. Empirically the r-ratio test will pin c: lim r(t)/t² = c, then expected err(1e-4) = c·1e-4. I'll compute c first, then test A-limit err ≈ c·t within 20%.\n\nOrder: compute r-ratio first to get c, then A-limit test against A + c·t.\n\nActually r(t)/t² → c, and g(t)/t − ½log t − A = r(t)/t = c·t + o(t). So test |observed − c·t| / |c·t| small. Good.\n\nBut note: the A-limit test uses t below ½ log 2 ≈ 0.3466; 1e-2..1e-4 fine. r-ratio uses t up to 0.25 (also below 0.3466 ✔).\n\nNow the table: use mpmath mp.dps = 25 for lerchphi, but prime sum from float cumsum (error ~1e-13·magnitude, fine at 25dps? It contaminates the 25-digit value — but we only print ~17 digits to CSV. Set dps 20 for speed.)\n\nRuntime estimate for 10001 lerchphi calls at dps 20: unknown; add time printout. If too slow I'll subsample near 0. Let me include a fallback: compute lerch combination via direct series when e^{−2t} < 0.5 (t > 0.35): fast numpy vectorized. For t ≤ 0.35 (350 points) use mpmath. Actually z = e^{−2t} < 0.5 means terms decay by ≥0.5^k — need only ~50 terms for 1e-15. For t>0.35 vectorized numpy over 10k points × 60 terms trivial. For t ≤ 0.35: 351 mpmath calls — fast. \n\nAlso at t=0: g=0 exactly, set directly.\n\nNumpy direct series: Phi(z) = Σ_{k=0}^{K} z^k/(k+0.25)^2. For array of z values (all < 0.5), K=80.\n\nThen combination C(t) = −0.25*(Phi1 − exp(−t/2)*Phi(exp(−2t))) where Phi1 = polygamma(1,0.25) float64 = 17.19732915450756... get via mpmath once at dps 30, hardcode as computed value in JSON.\n\nAlso −4(cosh term): −4*(exp(t/2)+exp(−t/2)−2).\n\ndigamma term float: coeff = 2.6860917096128... compute via mpmath once.\n\nPrime sum via cumsum: logs = log(n) for n=2..N; sieve Λ; lam = Λ/√n; C1 = cumsum(lam), C0 = cumsum(lam*log n); m = np.searchsorted(logs, t, 'right'); P = C1[m−1]*t − C0[m−1] (if m>0).\n\nFor high-precision validation prime sum: direct mpmath loop over prime powers ≤ e^t (t ≤ 0.25 → n ≤ 1, wait e^0.25 = 1.28 < 2, prime sum = 0!). Oh nice — for ALL validation t ≤ 0.3466 the prime sum vanishes. So validation needs no prime sum at all. But evenness test at t=2.3, 5.1 does need it — use mpmath loop there (n ≤ e^5.1 = 164, small).\n\nSieve: N = int(exp(10))+2 = 22028. lam array.\n\nNow write the script.\n\ndef von_mangoldt_sieve(N): returns array lam_log[0..N] where lam_log[n] = log p if n = p^k.\n```python\nlam_log = np.zeros(N+1)\nis_composite = np.zeros(N+1, bool)\nfor p in range(2, N+1):\n if not is_composite[p]:\n pk = p\n while pk <= N:\n lam_log[pk] = np.log(p)\n # mark multiples of p composite? standard sieve: mark p-multiples\n pk *= p\n for m in range(2*p, N+1, p):\n is_composite[m] = True\n```\nWait marking must happen for all p; do standard sieve first for primes, then powers.\n\nCSV rows 10001 — fine.\n\nJSON checks and key values.\n\nLet me also add g evaluation at interesting points for the record: g(½log2) (first prime kink), g(log 2), g(log 4)... record in JSON.\n\nWrite the file now.","signature":""},{"type":"text","text":"Conventions clear. Writing the E-009a script now."},{"type":"tool_use","id":"Write_8","name":"Write","input":{"file_path":"D:\\Python\\Quantum\\research\\numerics\\scripts\\E-009-suzuki-screw-function.py","content":"\"\"\"E-009a — Task 3.4c (Q-3.4c, D-018 plan step 1): Suzuki's screw function g(t).\n\nSuzuki (arXiv:2606.09096v1, eq. 1.3) — explicit screw function of zeta:\n\n g(t) = -4(e^{t/2} + e^{-t/2} - 2)\n + Sum_{n <= e^{|t|}} (Lambda(n)/sqrt(n)) (|t| - log n)\n - (|t|/2)(psi(1/4) - log pi)\n - (1/4)[ Phi(1,2,1/4) - e^{-|t|/2} Phi(e^{-2|t|}, 2, 1/4) ]\n\n(psi = digamma, Phi = Hurwitz-Lerch transcendent Phi(z,s,a) =\nSum_{k>=0} z^k/(k+a)^s; Lambda = von Mangoldt; prime sum = 0 for\n|t| < (1/2) log 2.)\n\nNear-origin expansion (Suzuki eq. 2.2):\n g(t) = (1/2)|t| log|t| + A|t| + (prime sum) + r(t),\n r even C^2, r = O(t^2),\n A = (1/2)(log(2pi) - psi(2)) = (1/2)(log(2pi) + gamma_E - 1)\n = 0.7075463656...\n\ng is a Krein-Langer screw function IFF RH (this is the Cor 1.6 bridge).\n\nPart A (validation, mpmath dps 30):\n A1. Lerch convention: Phi(1,2,1/4) == psi'(1/4) = polygamma(1, 1/4).\n A2. Evenness: g(-t) == g(t) at t = 0.7, 2.3, 5.1 (exercises prime sum).\n A3. A-limit: g(t)/t - (1/2) log t -> A as t -> 0+ (prime sum = 0 there);\n residual must equal r(t)/t = c2*t + o(t) with c2 from A4.\n A4. r-curvature: r(t)/t^2 -> c2 (constant) as t -> 0+, t <= 0.25.\n\nPart B (tabulation, float64 + mpmath hybrid):\n g(t) on t in [0, 10], step 1e-3 (mesh for the E-009b Galerkin stage;\n g is even so t >= 0 suffices). Prime-sum part via exact cumsum arrays\n (P(t) piecewise linear in t with kinks at log n). Lerch combination:\n direct numpy series for t > 0.35 (z = e^{-2t} < 0.5, fast), mpmath\n lerchphi for t <= 0.35 (slow-converging near z = 1).\n\nOutputs: research/numerics/data/E-009-suzuki-screw-function.csv (+ .json)\n\"\"\"\nimport json\nimport sys\nimport time\nfrom datetime import datetime, timezone\nfrom pathlib import Path\n\nimport numpy as np\nimport mpmath as mp\n\nREPO = Path(__file__).resolve().parents[3]\nDATA = REPO / \"research\" / \"numerics\" / \"data\"\nCSV = DATA / \"E-009-suzuki-screw-function.csv\"\nJSON = DATA / \"E-009-suzuki-screw-function.json\"\n\nT_MAX = 10.0 # tabulation range [0, T_MAX] (g even)\nH = 1e-3 # mesh step\nDPS = 30 # validation precision\nA_CONST = 0.5 * (np.log(2 * np.pi) + np.euler_gamma - 1.0) # eq. 2.2 constant\n\n\n# ---------------------------------------------------------------- von Mangoldt\ndef mangoldt_sieve(n_max: int) -> np.ndarray:\n \"\"\"lam_log[n] = log p if n = p^k (prime power), else 0.\"\"\"\n lam_log = np.zeros(n_max + 1)\n is_comp = np.zeros(n_max + 1, dtype=bool)\n for p in range(2, n_max + 1):\n if not is_comp[p]:\n if p * p <= n_max:\n is_comp[p * p :: p] = True\n pk = p\n while pk <= n_max:\n lam_log[pk] = np.log(p)\n pk *= p\n return lam_log\n\n\nN_MAX = int(np.exp(T_MAX)) + 2\nLAM_LOG = mangoldt_sieve(N_MAX)\nNS = np.arange(2, N_MAX + 1)\nLOGS = np.log(NS)\nLAM = LAM_LOG[2:] / np.sqrt(NS) # Lambda(n)/sqrt(n)\nC1 = np.cumsum(LAM) # Sum_{n<=x} Lambda(n)/sqrt(n)\nC0 = np.cumsum(LAM * LOGS) # Sum_{n<=x} Lambda(n) log(n)/sqrt(n)\n\n\ndef prime_sum_vec(t: np.ndarray) -> np.ndarray:\n \"\"\"P(t) = Sum_{n <= e^t} (Lambda(n)/sqrt(n))(t - log n), vectorized.\"\"\"\n m = np.searchsorted(LOGS, t, side=\"right\") # count of n with log n <= t\n out = np.zeros_like(t)\n pos = m > 0\n idx = m[pos] - 1\n out[pos] = C1[idx] * t[pos] - C0[idx]\n return out\n\n\ndef prime_sum_mp(t: mp.mpf) -> mp.mpf:\n s = mp.mpf(0)\n n = 2\n while mp.log(n) <= t:\n s += (LAM_LOG[n] / mp.sqrt(n)) * (t - mp.log(n))\n n += 1\n return s\n\n\n# ------------------------------------------------------------------- g(t), mp\ndef g_mp(t) -> mp.mpf:\n t = mp.mpf(t)\n a = abs(t)\n s = -4 * (mp.exp(a / 2) + mp.exp(-a / 2) - 2)\n s += prime_sum_mp(a)\n s -= (a / 2) * (mp.digamma(mp.mpf(1) / 4) - mp.log(mp.pi))\n s -= mp.mpf(1) / 4 * (mp.lerchphi(1, 2, mp.mpf(1) / 4)\n - mp.exp(-a / 2) * mp.lerchphi(mp.exp(-2 * a), 2,\n mp.mpf(1) / 4))\n return s\n\n\n# --------------------------------------------------------------- g(t), float\nKK = np.arange(0, 80)\nDENOM = (KK + 0.25) ** 2\n\n\ndef lerch_phi2_fast(z: np.ndarray) -> np.ndarray:\n \"\"\"Phi(z,2,1/4) by direct series; valid for |z| <= 0.5 (t >= 0.35).\"\"\"\n powers = np.power.outer(z, KK) # (n_z, K)\n return (powers / DENOM).sum(axis=1)\n\n\ndef g_float_table(t: np.ndarray) -> np.ndarray:\n g = -4.0 * (np.exp(t / 2) + np.exp(-t / 2) - 2.0)\n g += prime_sum_vec(t)\n dig_coeff = -0.5 * float(mp.digamma(mp.mpf(1) / 4) - mp.log(mp.pi))\n g += dig_coeff * t\n phi1 = float(mp.lerchphi(1, 2, mp.mpf(1) / 4)) # = polygamma(1, 1/4)\n lerch_comb = np.empty_like(t)\n slow = t <= 0.35\n if slow.any():\n mp.dps = 25\n for i in np.flatnonzero(slow):\n ti = mp.mpf(float(t[i]))\n val = (mp.lerchphi(1, 2, mp.mpf(1) / 4)\n - mp.exp(-ti / 2) * mp.lerchphi(mp.exp(-2 * ti), 2,\n mp.mpf(1) / 4))\n lerch_comb[i] = float(val) / -4.0\n mp.dps = DPS\n if (~slow).any():\n z = np.exp(-2.0 * t[~slow])\n lerch_comb[~slow] = -0.25 * (phi1 - np.exp(-t[~slow] / 2)\n * lerch_phi2_fast(z))\n g += lerch_comb\n g[t == 0.0] = 0.0\n return g\n\n\n# ------------------------------------------------------------------- main\ndef main() -> int:\n if CSV.exists() and \"--force\" not in sys.argv:\n print(f\"cached: {CSV} exists (use --force to recompute)\"); return 0\n t0 = time.time()\n mp.dps = DPS\n\n # ---- A1: Lerch convention\n phi1_mp = mp.lerchphi(1, 2, mp.mpf(1) / 4)\n trig1_mp = mp.polygamma(1, mp.mpf(1) / 4)\n a1_err = abs(phi1_mp - trig1_mp)\n check_a1 = bool(a1_err < mp.mpf(\"1e-25\"))\n\n # ---- A2: evenness (t values exercising the prime sum)\n ev_ts = [mp.mpf(\"0.7\"), mp.mpf(\"2.3\"), mp.mpf(\"5.1\")]\n ev_err = max(abs(g_mp(t) - g_mp(-t)) for t in ev_ts)\n check_a2 = bool(ev_err < mp.mpf(\"1e-20\"))\n\n # ---- A4: r(t)/t^2 -> c2 (t <= 0.25: prime sum = 0)\n A_mp = (mp.log(2 * mp.pi) + mp.euler - 1) / 2\n r_ts = [mp.mpf(\"0.25\"), mp.mpf(\"0.125\"), mp.mpf(\"0.0625\"),\n mp.mpf(\"0.03125\"), mp.mpf(\"0.015625\")]\n r_ratios = []\n for t in r_ts:\n r = g_mp(t) - (t / 2) * mp.log(t) - A_mp * t # prime sum = 0 here\n r_ratios.append(r / t**2)\n c2 = r_ratios[-1]\n r_gap_last = abs(r_ratios[-1] - r_ratios[-2])\n r_gap_prev = abs(r_ratios[-2] - r_ratios[-3])\n check_a4 = bool(r_gap_last < r_gap_prev) # gaps shrinking => converging\n\n # ---- A3: g(t)/t - (1/2) log t -> A; residual ~ c2 * t\n a3_ts = [mp.mpf(\"1e-2\"), mp.mpf(\"1e-3\"), mp.mpf(\"1e-4\")]\n a3_resid = [g_mp(t) / t - mp.log(t) / 2 - A_mp for t in a3_ts]\n # expected residual ~ c2*t; test relative agreement within 25%\n a3_rel = [abs((res - c2 * t) / (c2 * t)) if c2 != 0 else abs(res)\n for res, t in zip(a3_resid, a3_ts)]\n check_a3 = bool(a3_rel[-1] < 0.25 and abs(a3_resid[-1]) < mp.mpf(\"3e-3\"))\n\n print(f\"A1 lerch==polygamma(1,1/4): err {mp.nstr(a1_err, 5)}\")\n print(f\"A2 evenness max err: {mp.nstr(ev_err, 5)}\")\n print(f\"A4 r(t)/t^2 -> c2 = {mp.nstr(c2, 12)} \"\n f\"(gaps {mp.nstr(r_gap_prev, 3)} -> {mp.nstr(r_gap_last, 3)})\")\n print(f\"A3 residuals g(t)/t - (1/2)log t - A at 1e-2/1e-3/1e-4: \"\n f\"{[mp.nstr(x, 6) for x in a3_resid]} (c2*t predicted)\")\n\n # ---- Part B: table on [0, T_MAX]\n t_mesh = np.arange(0.0, T_MAX + 0.5 * H, H)\n g_vals = g_float_table(t_mesh)\n p_vals = prime_sum_vec(t_mesh)\n check_b = bool(np.all(np.isfinite(g_vals)))\n # spot check: float table vs mpmath at a few points\n spot_ts = [0.123, 1.0, 3.7, 7.0, 10.0]\n spot_err = max(abs(g_mp(x) - mp.mpf(float(g_vals[int(round(x / H))])))\n for x in spot_ts)\n check_spot = bool(spot_err < mp.mpf(\"1e-9\"))\n print(f\"Part B: {len(t_mesh)} points, spot-check err vs mpmath: \"\n f\"{mp.nstr(spot_err, 3)}\")\n\n checks = {\n \"A1_lerch_convention\": check_a1,\n \"A2_evenness\": check_a2,\n \"A3_A_limit\": check_a3,\n \"A4_r_curvature_converges\": check_a4,\n \"B_table_finite\": check_b,\n \"B_spot_check_vs_mpmath\": check_spot,\n }\n\n DATA.mkdir(parents=True, exist_ok=True)\n with CSV.open(\"w\", encoding=\"utf-8\", newline=\"\") as f:\n f.write(\"t,g_t,prime_sum_t\\n\")\n for ti, gi, pi in zip(t_mesh, g_vals, p_vals):\n f.write(f\"{ti:.4f},{gi:.12e},{pi:.12e}\\n\")\n\n meta = {\n \"experiment\": \"E-009a\",\n \"task\": \"3.4c (D-018 step 1): Suzuki screw function g(t)\",\n \"description\": \"mpmath implementation of Suzuki's explicit screw \"\n \"function g(t) (arXiv:2606.09096v1 eq. 1.3); \"\n \"validated against eq. 2.2 near-origin expansion; \"\n \"tabulated on [0,10] step 1e-3 for E-009b Galerkin\",\n \"script\": \"research/numerics/scripts/E-009-suzuki-screw-function.py\",\n \"parameters\": {\"t_max\": T_MAX, \"step\": H, \"validation_dps\": DPS,\n \"n_sieve_max\": N_MAX},\n \"checks\": checks,\n \"key_values\": {\n \"A_eq2.2\": float(A_CONST),\n \"c2_lim_r_over_t2\": float(c2),\n \"digamma_coeff_-half(psi(1/4)-log pi)\":\n float(-(mp.digamma(mp.mpf(1) / 4) - mp.log(mp.pi)) / 2),\n \"Phi(1,2,1/4)\": float(phi1_mp),\n \"g(log2)\": float(g_mp(mp.log(2))),\n \"g(1)\": float(g_mp(1)),\n \"g(5)\": float(g_mp(5)),\n \"g(10)\": float(g_vals[-1]),\n },\n \"caveats\": [\n \"float64 table: Lerch combination has cancellation for small t \"\n \"(~4 digits lost at t=1e-3); spot-checked vs mpmath to 1e-9\",\n \"prime sum is piecewise linear in t (kinks at log n) — exact \"\n \"via cumsum, no interpolation error\",\n \"g even: only t >= 0 tabulated\",\n ],\n \"elapsed_seconds\": round(time.time() - t0, 1),\n \"timestamp_utc\": datetime.now(timezone.utc).isoformat(),\n \"row_count\": int(len(t_mesh)),\n }\n JSON.write_text(json.dumps(meta, indent=2), encoding=\"utf-8\")\n ok = all(checks.values())\n print(f\"wrote {CSV.name}, {JSON.name} — \"\n f\"{'ALL CHECKS PASS' if ok else 'CHECK FAILURE: ' + str(checks)}\")\n return 0 if ok else 1\n\n\nif __name__ == \"__main__\":\n raise SystemExit(main())\n"},"extra_content":null}]