← Session 7 — All-forks plan approved; E-016 Herglotz work
thinking (600 chars)
E-014 ran in only 20s (sieve to e^20 apparently fast with the vectorized sieve... wait, sieve to e^20 = 4.85e8 — hmm, 20s total? The sieve must be fast with the vectorized approach, or... Actually N_MAXG = e^20 ≈ 4.85e8, arrays ~4GB each; on this machine maybe 60-90s. Whatever — 20.4s total says it's fast).
Hmm wait, actually T_MAXG=20 sieve might dominate. 20s total — fine. E-016 will be comparable. Foreground run OK (under 2-min default limit? 20s + derivative work — probably ~1 min; but E-015 is hogging 14 cores. Might take 3-5 min. Set timeout 600000 to be safe).
Now write E-016 script.
E-014 ran in 20s, so E-016 can run in the foreground. Writing the script now.
tool: Write ok
input
{
"file_path": "D:\\Python\\Quantum\\research\\numerics\\scripts\\E-016-herglotz-fine-structure.py",
"content": "\"\"\"E-016 \u2014 Q-S3c: Herglotz (Poisson) fine structure of W at theta=pi vs F.\n\nTheory (D-023 \u00a74, item 3): for entire W with ALL REAL zeros z_j,\n -Im(W'/W)(iy) = -Im(kappa) + sum_j y/(y^2 + z_j^2)\nis (pi times) the Poisson-smoothed zero density at height y. Individual\nzeros within |z_j| <~ y leave distinguishable fingerprints, so this\ndiagnostic detects sub-lattice-scale zero redistribution BEFORE the\ndensity collapses (Cor 1.6 requires the a->inf zeros to become the\ngamma_j / xi'-zeros of F = z^2 xi(1/2-iz)/xi'(1/2-iz)).\n\nFor F on the imaginary axis (z = iy, s = 1/2+y real):\n rho_F(y) := -(1/pi) Im F'/F(iy)\n = (1/pi)[ 2/y + xi'/xi(1/2+y) - xi''/xi'(1/2+y) ].\nEquivalently rho_F(y) = (1/pi) d/dy log|F(iy)| (finite-diff check V2).\n\nDiagnostics per a in {4,6,8,10} (h=0.01, lambda=-1, theta=pi; theta=0\ncontrol):\n D1. rho_W(y; a, pi) on y in [0.5, 20] via the ANALYTIC z-derivative of\n the Galerkin Fourier representation (no finite differences):\n V(z) = hat(z) sum c_j e^{i z x_j}\n V'(z) = hat'(z) sum c_j e^{i z x_j} + hat(z) sum (i x_j) c_j e^{i z x_j}\n W(z) = (z-i) Vp + e^{i theta} (z+i) Vm (and z-derivative).\n D2. rho_F(y) analytic; cross-checked vs finite differences of\n log|F(iy)| (V2).\n D3. Zero-Poisson reconstruction: zeros of W(.,pi) located on the real\n axis (|W|-minima polished); rho from the Poisson sum must match D1\n up to the additive constant -Im(kappa)/pi (mean-matched over the\n window) \u2014 validates the derivative machinery (V3).\n D4. Metrics vs a: plateau value vs a/pi (density law); L2 distance\n d(a,theta) = ||rho_W - rho_F|| on [1,20] (theta=pi vs theta=0 \u2014\n does theta=pi actually sharpen the approach?); oscillation\n envelope of rho_W at small y.\n\nChecks:\n V1. plateau: |mean rho_W(y in [10,20]) - a/pi| / (a/pi) < 0.05.\n V2. rho_F analytic vs finite-diff of log|F|: max rel dev < 1e-6.\n V3. Poisson reconstruction from located zeros matches D1 (mean-matched\n L2 / rms < 0.02).\n V4. mirror symmetry v-(x) = v+(-x) rel err < 1e-6 (machinery regression).\n\nOutputs: research/numerics/data/E-016-herglotz-fine-structure.csv (+ .json)\n\"\"\"\nimport json\nimport sys\nimport time\nfrom datetime import datetime, timezone\nfrom pathlib import Path\n\nimport mpmath as mp\nimport numpy as np\nimport scipy.linalg\nimport scipy.optimize\nfrom numpy.polynomial.legendre import leggauss\n\nREPO = Path(__file__).resolve().parents[3]\nDATA = REPO / \"research\" / \"numerics\" / \"data\"\nCSV = DATA / \"E-016-herglotz-fine-structure.csv\"\nJSON = DATA / \"E-016-herglotz-fine-structure.json\"\n\nT_MAXG = 20.0 # sieve to e^20 covers 2a <= 20 (a <= 10)\nT_LERCH = 0.05\nGL_N = 8\nA_SWEEP = [4.0, 6.0, 8.0, 10.0]\nH_MAIN = 0.01\nLAM_SHIFT = -1.0\nY_GRID = np.linspace(0.5, 20.0, 391)\nY_PLATEAU = (10.0, 20.0)\nZ_MAX_ZERO = 60.0\nN_GRID = 40000\n\n\n# ---------------------------------------------------------- g (as E-014)\ndef mangoldt_sieve_fast(n_max: int) -> np.ndarray:\n is_comp = np.zeros(n_max + 1, dtype=bool)\n for p in range(2, int(n_max ** 0.5) + 1):\n if not is_comp[p]:\n is_comp[p * p::p] = True\n lam_log = np.zeros(n_max + 1)\n primes = np.flatnonzero(~is_comp[2:]) + 2\n lp = np.log(primes)\n pk = primes.copy()\n while pk.size and pk[0] <= n_max:\n lam_log[pk] = lp\n pk = pk * primes[: pk.size]\n keep = pk <= n_max\n pk = pk[keep]\n lp = lp[keep]\n return lam_log\n\n\nprint(\"sieving von Mangoldt to e^20 ...\", flush=True)\nt_sv = time.time()\nN_MAXG = int(np.exp(T_MAXG)) + 2\nLAM_LOG = mangoldt_sieve_fast(N_MAXG)\nNS = np.arange(2, N_MAXG + 1, dtype=np.float64)\nLAM = LAM_LOG[2:] / np.sqrt(NS)\nLOG_NS = np.log(NS)\nC1 = np.cumsum(LAM)\nC0 = np.cumsum(LAM * LOG_NS)\ndel LAM, LOG_NS, NS\nprint(f\" sieve done: {time.time()-t_sv:.0f}s\", flush=True)\n\nmp.mp.dps = 30\nDIG_COEFF = float(-(mp.digamma(mp.mpf(1) / 4) - mp.log(mp.pi)) / 2)\nPHI1 = float(mp.lerchphi(1, 2, mp.mpf(1) / 4))\nKK = np.arange(0, 400)\nDENOM = (KK + 0.25) ** 2\n\n\ndef prime_sum_vec(t):\n m = np.floor(np.exp(t) + 1e-9).astype(np.int64) - 1\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 lerch_comb(t):\n out = np.empty_like(t)\n fast = t > T_LERCH\n if fast.any():\n z = np.exp(-2.0 * t[fast])\n phi = (np.power.outer(z, KK) / DENOM).sum(axis=1)\n out[fast] = -0.25 * (PHI1 - np.exp(-t[fast] / 2) * phi)\n if (~fast).any():\n for i in np.flatnonzero(~fast):\n ti = mp.mpf(float(t[i]))\n out[i] = float(-0.25 * (mp.lerchphi(1, 2, mp.mpf(1) / 4)\n - mp.exp(-ti / 2)\n * mp.lerchphi(mp.exp(-2 * ti), 2,\n mp.mpf(1) / 4)))\n return out\n\n\ndef g_float(t):\n t = np.asarray(t, dtype=float)\n g = -4.0 * (np.exp(t / 2) + np.exp(-t / 2) - 2.0)\n g += prime_sum_vec(t)\n g += DIG_COEFF * t\n g += lerch_comb(t)\n g[t == 0.0] = 0.0\n return g\n\n\nA_CONST = 0.5 * (np.log(2 * np.pi) + np.euler_gamma - 1.0)\n\n\ndef sigma_vec(t):\n t = np.asarray(t, dtype=float)\n with np.errstate(divide=\"ignore\", invalid=\"ignore\"):\n return np.where(t > 0, 0.5 * t * np.log(t) + A_CONST * t, 0.0)\n\n\ndef F_sing(t):\n if t <= 0:\n return 0.0\n return (A_CONST * t**3 / 6.0 + t**3 * np.log(t) / 12.0\n - 5.0 * t**3 / 72.0)\n\n\ndef S_sing(k, h):\n return (F_sing((k + 1) * h) - 2.0 * F_sing(k * h)\n + F_sing(abs(k - 1) * h))\n\n\nXN, XW = leggauss(GL_N)\n\n\ndef S_full(k, h):\n def quad_plain(lo, hi, peak):\n tot = 0.0\n for a, b in ((lo, peak), (peak, hi)):\n if b - a <= 0:\n continue\n mid, half = 0.5 * (a + b), 0.5 * (b - a)\n xs = mid + half * XN\n rho = g_float(xs) - sigma_vec(xs)\n tot += half * float(np.dot(XW, (h - np.abs(xs - peak)) * rho))\n return tot\n\n def quad(lo, hi, peak):\n n_lo = int(np.ceil(np.exp(lo)))\n n_hi = int(np.floor(np.exp(hi)))\n if n_hi - n_lo > 64:\n return quad_plain(lo, hi, peak)\n pts = sorted({lo, peak, hi} |\n {float(np.log(n))\n for n in range(max(n_lo, 2), n_hi + 1)\n if lo < np.log(n) < hi})\n tot = 0.0\n for p0, p1 in zip(pts, pts[1:]):\n if p1 - p0 <= 0:\n continue\n mid, half = 0.5 * (p0 + p1), 0.5 * (p1 - p0)\n xs = mid + half * XN\n rho = g_float(xs) - sigma_vec(xs)\n tot += half * float(np.dot(XW, (h - np.abs(xs - peak)) * rho))\n return tot\n\n if k == 0:\n return S_sing(0, h) + 2.0 * quad(0.0, h, 0.0)\n return S_sing(k, h) + quad((k - 1) * h, (k + 1) * h, k * h)\n\n\ndef stiffness_first_col(n_elem, h):\n n_q = n_elem - 1\n s = np.array([S_full(k, h) for k in range(n_elem)])\n q = np.empty(n_q)\n q[0] = (2.0 * s[0] - 2.0 * s[1]) / h**2\n for k in range(1, n_q):\n q[k] = (2.0 * s[k] - s[k + 1] - s[k - 1]) / h**2\n return q\n\n\n_BUILD = {}\n\n\ndef build(a, h):\n if (a, h) in _BUILD:\n return _BUILD[(a, h)]\n n_elem = int(round(2 * a / h))\n n = n_elem - 1\n q = stiffness_first_col(n_elem, h)\n Q = scipy.linalg.toeplitz(q)\n M = scipy.linalg.toeplitz(np.r_[2 * h / 3, h / 6,\n np.zeros(max(n - 2, 0))])\n ii = np.arange(1, n_elem)\n x = -a + ii * h\n xg, wg = leggauss(8)\n b = {}\n for sgn in (+1.0, -1.0):\n mid1 = x - h / 2\n xs1 = mid1[:, None] + (h / 2) * xg[None, :]\n phi1 = (xs1 - (x - h)[:, None]) / h\n mid2 = x + h / 2\n xs2 = mid2[:, None] + (h / 2) * xg[None, :]\n phi2 = ((x + h)[:, None] - xs2) / h\n b[sgn] = ((h / 2) * (np.exp(sgn * xs1) * phi1) @ wg\n + (h / 2) * (np.exp(sgn * xs2) * phi2) @ wg)\n _BUILD[(a, h)] = (x, h, Q, M, b)\n return _BUILD[(a, h)]\n\n\ndef solve_v(Q, M, b_s, lam, n):\n T = Q - lam * M\n e = np.ones(n)\n B = np.zeros((n + 1, n + 1))\n B[:n, :n] = T\n B[:n, n] = e\n B[n, :n] = e\n sol = scipy.linalg.solve(B, np.r_[b_s, 0.0], assume_a=\"sym\")\n return sol[:n], T\n\n\n# ------------------------------------------- V(z), V'(z), W(z), W'(z)\ndef hat_FT(z, h):\n return h * np.sinc(z * h / (2 * np.pi)) ** 2\n\n\ndef hat_FT_d(z, h):\n \"\"\"d/dz [h sinc(zh/2pi)^2], sinc(x)=sin(pi x)/(pi x).\"\"\"\n u = z * h / 2.0\n s = np.sinc(z * h / (2 * np.pi)) # sin(u)/u\n with np.errstate(divide=\"ignore\", invalid=\"ignore\"):\n ds = np.where(np.abs(u) > 1e-8,\n (u * np.cos(u) - np.sin(u)) / u**2,\n -u / 3.0)\n return 2.0 * s * ds * (h / 2.0)\n\n\ndef V_and_dV(z, c, x, h):\n E = np.exp(1j * np.outer(z, x))\n S1 = E @ c\n S2 = (1j * x * E) @ c\n return hat_FT(z, h) * S1, hat_FT_d(z, h) * S1 + hat_FT(z, h) * S2\n\n\ndef W_and_dW(z, Vp, dVp, Vm, dVm, theta):\n eit = np.exp(1j * theta)\n W = (z - 1j) * Vp + eit * (z + 1j) * Vm\n dW = Vp + (z - 1j) * dVp + eit * (Vm + (z + 1j) * dVm)\n return W, dW\n\n\n# --------------------------------------------------------------- F side\ndef xi_real(s):\n s = mp.mpf(s)\n return 0.5 * s * (s - 1) * mp.pi ** (-s / 2) * mp.gamma(s / 2) * mp.zeta(s)\n\n\nmp.mp.dps = 40\n\n\ndef rho_F_analytic(y):\n \"\"\"rho_F(y) = (1/pi)[2/y + xi'/xi(1/2+y) - xi''/xi'(1/2+y)].\"\"\"\n s = mp.mpf(0.5) + mp.mpf(y)\n f = lambda t: xi_real(t)\n x0 = xi_real(s)\n x1 = mp.diff(f, s)\n x2 = mp.diff(f, s, 2)\n return float((2.0 / mp.mpf(y) + x1 / x0 - x2 / x1) / mp.pi)\n\n\ndef logabs_F_imag(y):\n s = mp.mpf(0.5) + mp.mpf(y)\n xiv = xi_real(s)\n dx = mp.diff(xi_real, s)\n return mp.log(abs(-mp.mpf(y) ** 2 * xiv / dx))\n\n\ndef main() -> int:\n if CSV.exists() and \"--force\" not in sys.argv:\n print(f\"cached: {CSV} (use --force to recompute)\")\n return 0\n t0 = time.time()\n\n # ---- F side\n print(\"F side: rho_F analytic + finite-diff cross-check ...\", flush=True)\n rhoF = np.array([rho_F_analytic(y) for y in Y_GRID])\n dlogF = np.array([\n float((logabs_F_imag(mp.mpf(y) + 1e-4)\n - logabs_F_imag(mp.mpf(y) - 1e-4)) / mp.mpf(2e-4))\n for y in [\"%.10f\" % y for y in Y_GRID]]) / np.pi\n rel = np.abs(rhoF - dlogF) / np.maximum(np.abs(rhoF), 1e-12)\n v2_err = float(rel.max())\n print(f\" V2 finite-diff max rel dev: {v2_err:.2e}\", flush=True)\n\n rows = []\n for y, r in zip(Y_GRID, rhoF):\n rows.append(f\"F,{y:.6f},{r:.10e}\")\n\n sym_err = 0.0\n metrics = {}\n for a in A_SWEEP:\n print(f\"a={a}: solve v\u00b1, Herglotz densities ...\", flush=True)\n x, hh, Q, M, b = build(a, H_MAIN)\n n = len(x)\n cp, T = solve_v(Q, M, b[+1.0], LAM_SHIFT, n)\n cm, _ = solve_v(Q, M, b[-1.0], LAM_SHIFT, n)\n np_ = float(cp @ T @ cp)\n nm = float(cm @ T @ cm)\n cm = cm * np.sqrt(np_ / nm)\n sym_err = max(sym_err,\n float(np.max(np.abs(cm - cp[::-1]))\n / max(1.0, float(np.max(np.abs(cp))))))\n\n ziy = 1j * Y_GRID\n Vp, dVp = V_and_dV(ziy, cp, x, hh)\n Vm, dVm = V_and_dV(ziy, cm, x, hh)\n\n out = {}\n for th, tag in ((np.pi, \"pi\"), (0.0, \"zero\")):\n W, dW = W_and_dW(ziy, Vp, dVp, Vm, dVm, th)\n rhoW = -np.imag(dW / W) / np.pi\n out[tag] = rhoW\n for y, r in zip(Y_GRID, rhoW):\n rows.append(f\"W_th{tag},{a},{y:.6f},{r:.10e}\")\n\n # ---- zeros of W(., pi) on the real axis (for V3)\n zg = np.linspace(1e-3, Z_MAX_ZERO, N_GRID)\n Vpr, _ = V_and_dV(zg, cp, x, hh)\n Vmr, _ = V_and_dV(zg, cm, x, hh)\n Wr, _ = W_and_dW(zg, Vpr, None, Vmr, None, np.pi)\n aw = np.abs(Wr)\n loc = np.flatnonzero((aw[1:-1] < aw[:-2]) & (aw[1:-1] < aw[2:])) + 1\n zeros = []\n for i in loc:\n try:\n r = scipy.optimize.minimize_scalar(\n lambda zz: abs(W_and_dW(\n np.array([zz]),\n *V_and_dV(np.array([zz]), cp, x, hh),\n *V_and_dV(np.array([zz]), cm, x, hh),\n np.pi)[0][0]),\n bracket=(zg[i - 1], zg[i], zg[i + 1]),\n method=\"brent\",\n options={\"xtol\": 1e-12})\n if r.success and abs(r.fun) < 1e-6 * max(1.0, aw[i]):\n zeros.append(r.x)\n except Exception:\n pass\n zeros = np.array(zeros)\n\n # V3: Poisson reconstruction (+- symmetric zeros), mean-matched\n recon = np.empty_like(Y_GRID)\n for k, y in enumerate(Y_GRID):\n recon[k] = (2.0 * np.sum(y / (y**2 + zeros**2))) / np.pi\n rhoW = out[\"pi\"]\n shift = float(np.mean(rhoW - recon))\n rms = float(np.sqrt(np.mean(rhoW**2)))\n v3_rel = float(np.sqrt(np.mean((rhoW - recon - shift)**2)) / rms)\n\n # D4 metrics\n pm = (Y_GRID >= Y_PLATEAU[0]) & (Y_GRID <= Y_PLATEAU[1])\n plateau = float(np.mean(rhoW[pm]))\n dens = a / np.pi\n sel = Y_GRID >= 1.0\n d_pi = float(np.sqrt(np.mean((out[\"pi\"][sel] - rhoF[sel])**2)))\n d_zero = float(np.sqrt(np.mean((out[\"zero\"][sel] - rhoF[sel])**2)))\n metrics[a] = {\n \"n_zeros\": len(zeros), \"z1\": float(zeros[0]) if len(zeros) else None,\n \"plateau\": plateau, \"density_a_over_pi\": dens,\n \"plateau_rel_dev\": abs(plateau - dens) / dens,\n \"v3_rel\": v3_rel, \"kappa_shift\": shift,\n \"d_pi\": d_pi, \"d_zero\": d_zero,\n }\n print(f\" zeros={len(zeros)} z1={metrics[a]['z1']:.4f} \"\n f\"plateau={plateau:.4f} (a/pi={dens:.4f}) \"\n f\"V3 rel={v3_rel:.2e} d_pi={d_pi:.4f} d_zero={d_zero:.4f}\",\n flush=True)\n\n checks = {\n \"V1_plateau_density\": bool(all(\n metrics[a][\"plateau_rel_dev\"] < 0.05 for a in A_SWEEP)),\n \"V2_rhoF_finitediff\": bool(v2_err < 1e-6),\n \"V3_poisson_reconstruction\": bool(all(\n metrics[a][\"v3_rel\"] < 0.02 for a in A_SWEEP)),\n \"V4_vmirror_symmetry\": bool(sym_err < 1e-6),\n }\n\n with CSV.open(\"w\", encoding=\"utf-8\", newline=\"\") as f:\n f.write(\"series,a,y,rho\\n\")\n for r in rows:\n parts = r.split(\",\")\n if parts[0] == \"F\":\n f.write(f\"F,-,{parts[1]},{parts[2]}\\n\")\n else:\n f.write(f\"{parts[0]},{parts[1]},{parts[2]},{parts[3]}\\n\")\n\n meta = {\n \"experiment\": \"E-016\",\n \"task\": \"Q-S3c: Herglotz fine structure of W at theta=pi vs F\",\n \"description\": \"Poisson-smoothed zero density -(1/pi) Im W'/W(iy) \"\n \"at theta=pi (and theta=0 control) vs rho_F; \"\n \"analytic Galerkin z-derivative; zero-Poisson \"\n \"reconstruction cross-check\",\n \"script\": \"research/numerics/scripts/E-016-herglotz-fine-structure.py\",\n \"parameters\": {\"a_sweep\": A_SWEEP, \"h\": H_MAIN, \"lambda\": LAM_SHIFT,\n \"y_grid\": [0.5, 20.0, 391], \"z_max_zeros\": Z_MAX_ZERO},\n \"checks\": checks,\n \"key_values\": {\n \"rhoF_finitediff_max_rel\": v2_err,\n \"metrics\": {str(a): metrics[a] for a in A_SWEEP},\n \"vmirror_sym_err\": sym_err,\n },\n \"caveats\": [\n \"Galerkin discretization h=0.01 (E-010 h-convergence: z1 rel \"\n \"drift 4e-4 between h=0.01 and 0.005 at a=4)\",\n \"Poisson smoothing at height y averages over scale ~y: lattice \"\n \"oscillations at spacing pi/a are suppressed ~e^{-2ay}; the \"\n \"sensitive region is y <~ few (low zeros resolve individually)\",\n \"theta=pi vs theta=0 comparison is the sharp test of the H4 \"\n \"read-off (D-023): if theta=pi is correct, d_pi should lie \"\n \"systematically below d_zero as a grows\",\n ],\n \"elapsed_seconds\": round(time.time() - t0, 1),\n \"timestamp_utc\": datetime.now(timezone.utc).isoformat(),\n \"row_count\": None,\n }\n meta[\"row_count\"] = sum(1 for _ in CSV.open())\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"
}
result 2026-08-18 10:44:15Z → 2026-08-27 16:15:57Z
{
"type": "create",
"filePath": "D:\\Python\\Quantum\\web\\dashboard\\templates\\dashboard\\search.html",
"content": "{% extends \"dashboard/base.html\" %}\n{% block title %}Search{% endblock %}\n{% block nav_search %}active{% endblock %}\n{% block content %}\n<h1>Search the research archive</h1>\n<p class=\"muted\">Case-insensitive search across all research markdown \u2014\nresults, derivations, graveyard, journal, resources atlas. Multiple words\nare ANDed (every word must appear in the file).</p>\n\n<form method=\"get\" class=\"search-form\">\n <input type=\"text\" name=\"q\" value=\"{{ q }}\" maxlength=\"200\" autofocus\n placeholder=\"e.g. dog-leg or partialFraction_identity\">\n <button type=\"submit\">Search</button>\n</form>\n\n{% if search %}\n {% if search.results %}\n <p class=\"muted\">{{ search.total }} file{{ search.total|pluralize }}\n matched\n {% if search.truncated %} \u2014 showing the first {{ search.total }}\n (narrow the query for more){% endif %}</p>\n {% for kind, hits in search.results.items %}\n <h2>{{ kind|capfirst }} <span class=\"muted\">({{ hits|length }})</span></h2>\n {% for h in hits %}\n <div class=\"card search-hit\">\n <div class=\"search-hit-head\">\n <a href=\"{{ h.url }}\">{{ h.name }}</a>\n {% if h.title != h.name %}<span class=\"muted\">\u2014 {{ h.title }}</span>{% endif %}\n </div>\n {% for line_no, before, match, after in h.snippets %}\n <div class=\"snippet\"><span class=\"snippet-ln\">L{{ line_no }}</span>\u2026{{ before }}<mark>{{ match }}</mark>{{ after }}\u2026</div>\n {% endfor %}\n </div>\n {% endfor %}\n {% endfor %}\n {% else %}\n <p>No files matched <b>{{ search.q }}</b>.</p>\n {% endif %}\n{% endif %}\n{% endblock %}\n",
"structuredPatch": [],
"originalFile": null
}
raw content_json (database row)
[{"type":"thinking","thinking":"E-014 ran in only 20s (sieve to e^20 apparently fast with the vectorized sieve... wait, sieve to e^20 = 4.85e8 — hmm, 20s total? The sieve must be fast with the vectorized approach, or... Actually N_MAXG = e^20 ≈ 4.85e8, arrays ~4GB each; on this machine maybe 60-90s. Whatever — 20.4s total says it's fast). \n\nHmm wait, actually T_MAXG=20 sieve might dominate. 20s total — fine. E-016 will be comparable. Foreground run OK (under 2-min default limit? 20s + derivative work — probably ~1 min; but E-015 is hogging 14 cores. Might take 3-5 min. Set timeout 600000 to be safe).\n\nNow write E-016 script.","signature":""},{"type":"text","text":"E-014 ran in 20s, so E-016 can run in the foreground. Writing the script now."},{"type":"tool_use","id":"Write_25","name":"Write","input":{"file_path":"D:\\Python\\Quantum\\research\\numerics\\scripts\\E-016-herglotz-fine-structure.py","content":"\"\"\"E-016 — Q-S3c: Herglotz (Poisson) fine structure of W at theta=pi vs F.\n\nTheory (D-023 §4, item 3): for entire W with ALL REAL zeros z_j,\n -Im(W'/W)(iy) = -Im(kappa) + sum_j y/(y^2 + z_j^2)\nis (pi times) the Poisson-smoothed zero density at height y. Individual\nzeros within |z_j| <~ y leave distinguishable fingerprints, so this\ndiagnostic detects sub-lattice-scale zero redistribution BEFORE the\ndensity collapses (Cor 1.6 requires the a->inf zeros to become the\ngamma_j / xi'-zeros of F = z^2 xi(1/2-iz)/xi'(1/2-iz)).\n\nFor F on the imaginary axis (z = iy, s = 1/2+y real):\n rho_F(y) := -(1/pi) Im F'/F(iy)\n = (1/pi)[ 2/y + xi'/xi(1/2+y) - xi''/xi'(1/2+y) ].\nEquivalently rho_F(y) = (1/pi) d/dy log|F(iy)| (finite-diff check V2).\n\nDiagnostics per a in {4,6,8,10} (h=0.01, lambda=-1, theta=pi; theta=0\ncontrol):\n D1. rho_W(y; a, pi) on y in [0.5, 20] via the ANALYTIC z-derivative of\n the Galerkin Fourier representation (no finite differences):\n V(z) = hat(z) sum c_j e^{i z x_j}\n V'(z) = hat'(z) sum c_j e^{i z x_j} + hat(z) sum (i x_j) c_j e^{i z x_j}\n W(z) = (z-i) Vp + e^{i theta} (z+i) Vm (and z-derivative).\n D2. rho_F(y) analytic; cross-checked vs finite differences of\n log|F(iy)| (V2).\n D3. Zero-Poisson reconstruction: zeros of W(.,pi) located on the real\n axis (|W|-minima polished); rho from the Poisson sum must match D1\n up to the additive constant -Im(kappa)/pi (mean-matched over the\n window) — validates the derivative machinery (V3).\n D4. Metrics vs a: plateau value vs a/pi (density law); L2 distance\n d(a,theta) = ||rho_W - rho_F|| on [1,20] (theta=pi vs theta=0 —\n does theta=pi actually sharpen the approach?); oscillation\n envelope of rho_W at small y.\n\nChecks:\n V1. plateau: |mean rho_W(y in [10,20]) - a/pi| / (a/pi) < 0.05.\n V2. rho_F analytic vs finite-diff of log|F|: max rel dev < 1e-6.\n V3. Poisson reconstruction from located zeros matches D1 (mean-matched\n L2 / rms < 0.02).\n V4. mirror symmetry v-(x) = v+(-x) rel err < 1e-6 (machinery regression).\n\nOutputs: research/numerics/data/E-016-herglotz-fine-structure.csv (+ .json)\n\"\"\"\nimport json\nimport sys\nimport time\nfrom datetime import datetime, timezone\nfrom pathlib import Path\n\nimport mpmath as mp\nimport numpy as np\nimport scipy.linalg\nimport scipy.optimize\nfrom numpy.polynomial.legendre import leggauss\n\nREPO = Path(__file__).resolve().parents[3]\nDATA = REPO / \"research\" / \"numerics\" / \"data\"\nCSV = DATA / \"E-016-herglotz-fine-structure.csv\"\nJSON = DATA / \"E-016-herglotz-fine-structure.json\"\n\nT_MAXG = 20.0 # sieve to e^20 covers 2a <= 20 (a <= 10)\nT_LERCH = 0.05\nGL_N = 8\nA_SWEEP = [4.0, 6.0, 8.0, 10.0]\nH_MAIN = 0.01\nLAM_SHIFT = -1.0\nY_GRID = np.linspace(0.5, 20.0, 391)\nY_PLATEAU = (10.0, 20.0)\nZ_MAX_ZERO = 60.0\nN_GRID = 40000\n\n\n# ---------------------------------------------------------- g (as E-014)\ndef mangoldt_sieve_fast(n_max: int) -> np.ndarray:\n is_comp = np.zeros(n_max + 1, dtype=bool)\n for p in range(2, int(n_max ** 0.5) + 1):\n if not is_comp[p]:\n is_comp[p * p::p] = True\n lam_log = np.zeros(n_max + 1)\n primes = np.flatnonzero(~is_comp[2:]) + 2\n lp = np.log(primes)\n pk = primes.copy()\n while pk.size and pk[0] <= n_max:\n lam_log[pk] = lp\n pk = pk * primes[: pk.size]\n keep = pk <= n_max\n pk = pk[keep]\n lp = lp[keep]\n return lam_log\n\n\nprint(\"sieving von Mangoldt to e^20 ...\", flush=True)\nt_sv = time.time()\nN_MAXG = int(np.exp(T_MAXG)) + 2\nLAM_LOG = mangoldt_sieve_fast(N_MAXG)\nNS = np.arange(2, N_MAXG + 1, dtype=np.float64)\nLAM = LAM_LOG[2:] / np.sqrt(NS)\nLOG_NS = np.log(NS)\nC1 = np.cumsum(LAM)\nC0 = np.cumsum(LAM * LOG_NS)\ndel LAM, LOG_NS, NS\nprint(f\" sieve done: {time.time()-t_sv:.0f}s\", flush=True)\n\nmp.mp.dps = 30\nDIG_COEFF = float(-(mp.digamma(mp.mpf(1) / 4) - mp.log(mp.pi)) / 2)\nPHI1 = float(mp.lerchphi(1, 2, mp.mpf(1) / 4))\nKK = np.arange(0, 400)\nDENOM = (KK + 0.25) ** 2\n\n\ndef prime_sum_vec(t):\n m = np.floor(np.exp(t) + 1e-9).astype(np.int64) - 1\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 lerch_comb(t):\n out = np.empty_like(t)\n fast = t > T_LERCH\n if fast.any():\n z = np.exp(-2.0 * t[fast])\n phi = (np.power.outer(z, KK) / DENOM).sum(axis=1)\n out[fast] = -0.25 * (PHI1 - np.exp(-t[fast] / 2) * phi)\n if (~fast).any():\n for i in np.flatnonzero(~fast):\n ti = mp.mpf(float(t[i]))\n out[i] = float(-0.25 * (mp.lerchphi(1, 2, mp.mpf(1) / 4)\n - mp.exp(-ti / 2)\n * mp.lerchphi(mp.exp(-2 * ti), 2,\n mp.mpf(1) / 4)))\n return out\n\n\ndef g_float(t):\n t = np.asarray(t, dtype=float)\n g = -4.0 * (np.exp(t / 2) + np.exp(-t / 2) - 2.0)\n g += prime_sum_vec(t)\n g += DIG_COEFF * t\n g += lerch_comb(t)\n g[t == 0.0] = 0.0\n return g\n\n\nA_CONST = 0.5 * (np.log(2 * np.pi) + np.euler_gamma - 1.0)\n\n\ndef sigma_vec(t):\n t = np.asarray(t, dtype=float)\n with np.errstate(divide=\"ignore\", invalid=\"ignore\"):\n return np.where(t > 0, 0.5 * t * np.log(t) + A_CONST * t, 0.0)\n\n\ndef F_sing(t):\n if t <= 0:\n return 0.0\n return (A_CONST * t**3 / 6.0 + t**3 * np.log(t) / 12.0\n - 5.0 * t**3 / 72.0)\n\n\ndef S_sing(k, h):\n return (F_sing((k + 1) * h) - 2.0 * F_sing(k * h)\n + F_sing(abs(k - 1) * h))\n\n\nXN, XW = leggauss(GL_N)\n\n\ndef S_full(k, h):\n def quad_plain(lo, hi, peak):\n tot = 0.0\n for a, b in ((lo, peak), (peak, hi)):\n if b - a <= 0:\n continue\n mid, half = 0.5 * (a + b), 0.5 * (b - a)\n xs = mid + half * XN\n rho = g_float(xs) - sigma_vec(xs)\n tot += half * float(np.dot(XW, (h - np.abs(xs - peak)) * rho))\n return tot\n\n def quad(lo, hi, peak):\n n_lo = int(np.ceil(np.exp(lo)))\n n_hi = int(np.floor(np.exp(hi)))\n if n_hi - n_lo > 64:\n return quad_plain(lo, hi, peak)\n pts = sorted({lo, peak, hi} |\n {float(np.log(n))\n for n in range(max(n_lo, 2), n_hi + 1)\n if lo < np.log(n) < hi})\n tot = 0.0\n for p0, p1 in zip(pts, pts[1:]):\n if p1 - p0 <= 0:\n continue\n mid, half = 0.5 * (p0 + p1), 0.5 * (p1 - p0)\n xs = mid + half * XN\n rho = g_float(xs) - sigma_vec(xs)\n tot += half * float(np.dot(XW, (h - np.abs(xs - peak)) * rho))\n return tot\n\n if k == 0:\n return S_sing(0, h) + 2.0 * quad(0.0, h, 0.0)\n return S_sing(k, h) + quad((k - 1) * h, (k + 1) * h, k * h)\n\n\ndef stiffness_first_col(n_elem, h):\n n_q = n_elem - 1\n s = np.array([S_full(k, h) for k in range(n_elem)])\n q = np.empty(n_q)\n q[0] = (2.0 * s[0] - 2.0 * s[1]) / h**2\n for k in range(1, n_q):\n q[k] = (2.0 * s[k] - s[k + 1] - s[k - 1]) / h**2\n return q\n\n\n_BUILD = {}\n\n\ndef build(a, h):\n if (a, h) in _BUILD:\n return _BUILD[(a, h)]\n n_elem = int(round(2 * a / h))\n n = n_elem - 1\n q = stiffness_first_col(n_elem, h)\n Q = scipy.linalg.toeplitz(q)\n M = scipy.linalg.toeplitz(np.r_[2 * h / 3, h / 6,\n np.zeros(max(n - 2, 0))])\n ii = np.arange(1, n_elem)\n x = -a + ii * h\n xg, wg = leggauss(8)\n b = {}\n for sgn in (+1.0, -1.0):\n mid1 = x - h / 2\n xs1 = mid1[:, None] + (h / 2) * xg[None, :]\n phi1 = (xs1 - (x - h)[:, None]) / h\n mid2 = x + h / 2\n xs2 = mid2[:, None] + (h / 2) * xg[None, :]\n phi2 = ((x + h)[:, None] - xs2) / h\n b[sgn] = ((h / 2) * (np.exp(sgn * xs1) * phi1) @ wg\n + (h / 2) * (np.exp(sgn * xs2) * phi2) @ wg)\n _BUILD[(a, h)] = (x, h, Q, M, b)\n return _BUILD[(a, h)]\n\n\ndef solve_v(Q, M, b_s, lam, n):\n T = Q - lam * M\n e = np.ones(n)\n B = np.zeros((n + 1, n + 1))\n B[:n, :n] = T\n B[:n, n] = e\n B[n, :n] = e\n sol = scipy.linalg.solve(B, np.r_[b_s, 0.0], assume_a=\"sym\")\n return sol[:n], T\n\n\n# ------------------------------------------- V(z), V'(z), W(z), W'(z)\ndef hat_FT(z, h):\n return h * np.sinc(z * h / (2 * np.pi)) ** 2\n\n\ndef hat_FT_d(z, h):\n \"\"\"d/dz [h sinc(zh/2pi)^2], sinc(x)=sin(pi x)/(pi x).\"\"\"\n u = z * h / 2.0\n s = np.sinc(z * h / (2 * np.pi)) # sin(u)/u\n with np.errstate(divide=\"ignore\", invalid=\"ignore\"):\n ds = np.where(np.abs(u) > 1e-8,\n (u * np.cos(u) - np.sin(u)) / u**2,\n -u / 3.0)\n return 2.0 * s * ds * (h / 2.0)\n\n\ndef V_and_dV(z, c, x, h):\n E = np.exp(1j * np.outer(z, x))\n S1 = E @ c\n S2 = (1j * x * E) @ c\n return hat_FT(z, h) * S1, hat_FT_d(z, h) * S1 + hat_FT(z, h) * S2\n\n\ndef W_and_dW(z, Vp, dVp, Vm, dVm, theta):\n eit = np.exp(1j * theta)\n W = (z - 1j) * Vp + eit * (z + 1j) * Vm\n dW = Vp + (z - 1j) * dVp + eit * (Vm + (z + 1j) * dVm)\n return W, dW\n\n\n# --------------------------------------------------------------- F side\ndef xi_real(s):\n s = mp.mpf(s)\n return 0.5 * s * (s - 1) * mp.pi ** (-s / 2) * mp.gamma(s / 2) * mp.zeta(s)\n\n\nmp.mp.dps = 40\n\n\ndef rho_F_analytic(y):\n \"\"\"rho_F(y) = (1/pi)[2/y + xi'/xi(1/2+y) - xi''/xi'(1/2+y)].\"\"\"\n s = mp.mpf(0.5) + mp.mpf(y)\n f = lambda t: xi_real(t)\n x0 = xi_real(s)\n x1 = mp.diff(f, s)\n x2 = mp.diff(f, s, 2)\n return float((2.0 / mp.mpf(y) + x1 / x0 - x2 / x1) / mp.pi)\n\n\ndef logabs_F_imag(y):\n s = mp.mpf(0.5) + mp.mpf(y)\n xiv = xi_real(s)\n dx = mp.diff(xi_real, s)\n return mp.log(abs(-mp.mpf(y) ** 2 * xiv / dx))\n\n\ndef main() -> int:\n if CSV.exists() and \"--force\" not in sys.argv:\n print(f\"cached: {CSV} (use --force to recompute)\")\n return 0\n t0 = time.time()\n\n # ---- F side\n print(\"F side: rho_F analytic + finite-diff cross-check ...\", flush=True)\n rhoF = np.array([rho_F_analytic(y) for y in Y_GRID])\n dlogF = np.array([\n float((logabs_F_imag(mp.mpf(y) + 1e-4)\n - logabs_F_imag(mp.mpf(y) - 1e-4)) / mp.mpf(2e-4))\n for y in [\"%.10f\" % y for y in Y_GRID]]) / np.pi\n rel = np.abs(rhoF - dlogF) / np.maximum(np.abs(rhoF), 1e-12)\n v2_err = float(rel.max())\n print(f\" V2 finite-diff max rel dev: {v2_err:.2e}\", flush=True)\n\n rows = []\n for y, r in zip(Y_GRID, rhoF):\n rows.append(f\"F,{y:.6f},{r:.10e}\")\n\n sym_err = 0.0\n metrics = {}\n for a in A_SWEEP:\n print(f\"a={a}: solve v±, Herglotz densities ...\", flush=True)\n x, hh, Q, M, b = build(a, H_MAIN)\n n = len(x)\n cp, T = solve_v(Q, M, b[+1.0], LAM_SHIFT, n)\n cm, _ = solve_v(Q, M, b[-1.0], LAM_SHIFT, n)\n np_ = float(cp @ T @ cp)\n nm = float(cm @ T @ cm)\n cm = cm * np.sqrt(np_ / nm)\n sym_err = max(sym_err,\n float(np.max(np.abs(cm - cp[::-1]))\n / max(1.0, float(np.max(np.abs(cp))))))\n\n ziy = 1j * Y_GRID\n Vp, dVp = V_and_dV(ziy, cp, x, hh)\n Vm, dVm = V_and_dV(ziy, cm, x, hh)\n\n out = {}\n for th, tag in ((np.pi, \"pi\"), (0.0, \"zero\")):\n W, dW = W_and_dW(ziy, Vp, dVp, Vm, dVm, th)\n rhoW = -np.imag(dW / W) / np.pi\n out[tag] = rhoW\n for y, r in zip(Y_GRID, rhoW):\n rows.append(f\"W_th{tag},{a},{y:.6f},{r:.10e}\")\n\n # ---- zeros of W(., pi) on the real axis (for V3)\n zg = np.linspace(1e-3, Z_MAX_ZERO, N_GRID)\n Vpr, _ = V_and_dV(zg, cp, x, hh)\n Vmr, _ = V_and_dV(zg, cm, x, hh)\n Wr, _ = W_and_dW(zg, Vpr, None, Vmr, None, np.pi)\n aw = np.abs(Wr)\n loc = np.flatnonzero((aw[1:-1] < aw[:-2]) & (aw[1:-1] < aw[2:])) + 1\n zeros = []\n for i in loc:\n try:\n r = scipy.optimize.minimize_scalar(\n lambda zz: abs(W_and_dW(\n np.array([zz]),\n *V_and_dV(np.array([zz]), cp, x, hh),\n *V_and_dV(np.array([zz]), cm, x, hh),\n np.pi)[0][0]),\n bracket=(zg[i - 1], zg[i], zg[i + 1]),\n method=\"brent\",\n options={\"xtol\": 1e-12})\n if r.success and abs(r.fun) < 1e-6 * max(1.0, aw[i]):\n zeros.append(r.x)\n except Exception:\n pass\n zeros = np.array(zeros)\n\n # V3: Poisson reconstruction (+- symmetric zeros), mean-matched\n recon = np.empty_like(Y_GRID)\n for k, y in enumerate(Y_GRID):\n recon[k] = (2.0 * np.sum(y / (y**2 + zeros**2))) / np.pi\n rhoW = out[\"pi\"]\n shift = float(np.mean(rhoW - recon))\n rms = float(np.sqrt(np.mean(rhoW**2)))\n v3_rel = float(np.sqrt(np.mean((rhoW - recon - shift)**2)) / rms)\n\n # D4 metrics\n pm = (Y_GRID >= Y_PLATEAU[0]) & (Y_GRID <= Y_PLATEAU[1])\n plateau = float(np.mean(rhoW[pm]))\n dens = a / np.pi\n sel = Y_GRID >= 1.0\n d_pi = float(np.sqrt(np.mean((out[\"pi\"][sel] - rhoF[sel])**2)))\n d_zero = float(np.sqrt(np.mean((out[\"zero\"][sel] - rhoF[sel])**2)))\n metrics[a] = {\n \"n_zeros\": len(zeros), \"z1\": float(zeros[0]) if len(zeros) else None,\n \"plateau\": plateau, \"density_a_over_pi\": dens,\n \"plateau_rel_dev\": abs(plateau - dens) / dens,\n \"v3_rel\": v3_rel, \"kappa_shift\": shift,\n \"d_pi\": d_pi, \"d_zero\": d_zero,\n }\n print(f\" zeros={len(zeros)} z1={metrics[a]['z1']:.4f} \"\n f\"plateau={plateau:.4f} (a/pi={dens:.4f}) \"\n f\"V3 rel={v3_rel:.2e} d_pi={d_pi:.4f} d_zero={d_zero:.4f}\",\n flush=True)\n\n checks = {\n \"V1_plateau_density\": bool(all(\n metrics[a][\"plateau_rel_dev\"] < 0.05 for a in A_SWEEP)),\n \"V2_rhoF_finitediff\": bool(v2_err < 1e-6),\n \"V3_poisson_reconstruction\": bool(all(\n metrics[a][\"v3_rel\"] < 0.02 for a in A_SWEEP)),\n \"V4_vmirror_symmetry\": bool(sym_err < 1e-6),\n }\n\n with CSV.open(\"w\", encoding=\"utf-8\", newline=\"\") as f:\n f.write(\"series,a,y,rho\\n\")\n for r in rows:\n parts = r.split(\",\")\n if parts[0] == \"F\":\n f.write(f\"F,-,{parts[1]},{parts[2]}\\n\")\n else:\n f.write(f\"{parts[0]},{parts[1]},{parts[2]},{parts[3]}\\n\")\n\n meta = {\n \"experiment\": \"E-016\",\n \"task\": \"Q-S3c: Herglotz fine structure of W at theta=pi vs F\",\n \"description\": \"Poisson-smoothed zero density -(1/pi) Im W'/W(iy) \"\n \"at theta=pi (and theta=0 control) vs rho_F; \"\n \"analytic Galerkin z-derivative; zero-Poisson \"\n \"reconstruction cross-check\",\n \"script\": \"research/numerics/scripts/E-016-herglotz-fine-structure.py\",\n \"parameters\": {\"a_sweep\": A_SWEEP, \"h\": H_MAIN, \"lambda\": LAM_SHIFT,\n \"y_grid\": [0.5, 20.0, 391], \"z_max_zeros\": Z_MAX_ZERO},\n \"checks\": checks,\n \"key_values\": {\n \"rhoF_finitediff_max_rel\": v2_err,\n \"metrics\": {str(a): metrics[a] for a in A_SWEEP},\n \"vmirror_sym_err\": sym_err,\n },\n \"caveats\": [\n \"Galerkin discretization h=0.01 (E-010 h-convergence: z1 rel \"\n \"drift 4e-4 between h=0.01 and 0.005 at a=4)\",\n \"Poisson smoothing at height y averages over scale ~y: lattice \"\n \"oscillations at spacing pi/a are suppressed ~e^{-2ay}; the \"\n \"sensitive region is y <~ few (low zeros resolve individually)\",\n \"theta=pi vs theta=0 comparison is the sharp test of the H4 \"\n \"read-off (D-023): if theta=pi is correct, d_pi should lie \"\n \"systematically below d_zero as a grows\",\n ],\n \"elapsed_seconds\": round(time.time() - t0, 1),\n \"timestamp_utc\": datetime.now(timezone.utc).isoformat(),\n \"row_count\": None,\n }\n meta[\"row_count\"] = sum(1 for _ in CSV.open())\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}]