E-001-window-constants.py

{
  "experiment": "E-001",
  "task": "0.4",
  "description": "Symbolic+numeric verification of R(psi0)=4/3 and R(psiMT)=1/2+(1/sqrt(2))cot(1/sqrt(2)); certificates 2-R and (3-R)/2",
  "script": "research/numerics/scripts/E-001-window-constants.py",
  "parameters": {
    "windows": [
      "psi0",
      "psiMT"
    ],
    "dps": 50
  },
  "symbolic_checks": {
    "R(psi0) == 4/3": true,
    "R(psiMT) == 1/2 + (1/sqrt2)cot(1/sqrt2) [60-digit]": true
  },
  "timestamp_utc": "2026-08-17T18:47:23.410351+00:00",
  "row_count": 2
}

"""E-001 — Task 0.4: symbolic + numeric verification of window constants.

R(psi) := [∫ψ² + ∬|u−v|ψ(u)ψ(v)dudv] / (∫ψ)²   on [−1/2, 1/2], ψ even.

Claims verified (paper Lemma 5.6; hand-derived in D-004 §7):
  R(ψ₀)   = 4/3                        (ψ₀ = indicator)
  R(ψ_MT) = 1/2 + (1/√2)cot(1/√2)      (ψ_MT = cos(√2 s))
and the certificates 2 − R(ψ), (3 − R(ψ))/2.

Outputs: research/numerics/data/E-001-window-constants.csv (+ .json sidecar)
"""
import json
import sys
from datetime import datetime, timezone
from pathlib import Path

import mpmath as mp
import sympy as sp

REPO = Path(__file__).resolve().parents[3]
DATA = REPO / "research" / "numerics" / "data"
DPS = 50

u, v = sp.symbols("u v", real=True)
half = sp.Rational(1, 2)


def R_symbolic(psi):
    """Compute R(ψ) symbolically for even ψ on [-1/2, 1/2]."""
    num1 = sp.integrate(psi**2, (u, -half, half))
    # ∬|u−v|ψ(u)ψ(v) = 2∫_{u>v}(u−v)ψ(u)ψ(v)  (integrand symmetric)
    inner = sp.integrate((u - v) * psi.subs(u, v), (v, -half, u))
    num2 = 2 * sp.integrate(sp.simplify(psi * inner), (u, -half, half))
    den = sp.integrate(psi, (u, -half, half)) ** 2
    return sp.simplify((num1 + num2) / den)


def main() -> int:
    mp.mp.dps = DPS
    ok = True

    psi0 = sp.Integer(1)
    psiMT = sp.cos(sp.sqrt(2) * u)

    R0 = R_symbolic(psi0)
    RMT = R_symbolic(psiMT)
    RMT_claim = sp.Rational(1, 2) + sp.cot(1 / sp.sqrt(2)) / sp.sqrt(2)

    checks = []
    checks.append(("R(psi0) == 4/3",
                   sp.simplify(R0 - sp.Rational(4, 3)) == 0, R0))
    # sympy's simplify does not fully reduce the trig form; verify the
    # difference vanishes to 50 digits instead of relying on simplify.
    diff = sp.N(RMT - RMT_claim, 60)
    checks.append(("R(psiMT) == 1/2 + (1/sqrt2)cot(1/sqrt2) [60-digit]",
                   abs(diff) < sp.Float("1e-50"), RMT))

    rows = []
    for name, R in [("psi0", sp.Rational(4, 3)), ("psiMT", RMT_claim)]:
        Rm = mp.mpf(sp.N(R, 60))
        rows.append({
            "window": name,
            "R_exact": str(R),
            "R_numeric": mp.nstr(Rm, 30),
            "2_minus_R": mp.nstr(2 - Rm, 30),
            "distinct_(3-R)/2": mp.nstr((3 - Rm) / 2, 30),
        })

    for label, passed, val in checks:
        print(f"[{'PASS' if passed else 'FAIL'}] {label}  (got {val})")
        ok &= passed

    print("\nCertificates:")
    for r in rows:
        print(f"  {r['window']:>6}: R = {r['R_numeric']}")
        print(f"         2 - R = {r['2_minus_R']}   (3-R)/2 = {r['distinct_(3-R)/2']}")

    csv_path = DATA / "E-001-window-constants.csv"
    json_path = DATA / "E-001-window-constants.json"
    DATA.mkdir(parents=True, exist_ok=True)
    if csv_path.exists() and "--force" not in sys.argv:
        print(f"\ncached: {csv_path} exists (use --force to recompute)")
        return 0 if ok else 1

    with csv_path.open("w", encoding="utf-8", newline="") as f:
        f.write(",".join(rows[0].keys()) + "\n")
        for r in rows:
            f.write(",".join(r.values()) + "\n")

    meta = {
        "experiment": "E-001",
        "task": "0.4",
        "description": "Symbolic+numeric verification of R(psi0)=4/3 and "
                       "R(psiMT)=1/2+(1/sqrt(2))cot(1/sqrt(2)); certificates "
                       "2-R and (3-R)/2",
        "script": "research/numerics/scripts/E-001-window-constants.py",
        "parameters": {"windows": ["psi0", "psiMT"], "dps": DPS},
        "symbolic_checks": {label: bool(passed) for label, passed, _ in checks},
        "timestamp_utc": datetime.now(timezone.utc).isoformat(),
        "row_count": len(rows),
    }
    json_path.write_text(json.dumps(meta, indent=2), encoding="utf-8")
    print(f"\nwrote {csv_path.name}, {json_path.name}")
    return 0 if ok else 1


if __name__ == "__main__":
    raise SystemExit(main())