E-013-xiprime-R-constants.py

{
  "experiment": "E-013",
  "task": "T4: xi'-analogue — why is the method stronger on xi'? (D-022)",
  "description": "R-functional with xi' diagonal law D1 (pair kernel x - 4x^2 + ...) reproduces Remark 7.1 constants; D1(x) <= x quantifies tight-pair suppression",
  "script": "research/numerics/scripts/E-013-xiprime-R-constants.py",
  "parameters": {
    "dps": 40,
    "d1_series_terms": 40
  },
  "checks": {
    "V1_flat_matches_085838": true,
    "V2_quartic_matches_086864": true,
    "V3_perfect_square_node_mechanism": true,
    "V4_xi_kernel_sanity_4_3": true
  },
  "key_values": {
    "R_xiprime_flat": "1.14161504167725",
    "cert_flat_computed": "0.858384958322751",
    "cert_flat_paper": 0.85838,
    "R_xiprime_quartic": "1.13135848298565",
    "cert_quartic_computed": "0.868641517014353",
    "cert_quartic_paper": 0.86864,
    "R_xi_flat_sanity": "1.33333242971717",
    "node_ratio_D1(0.5)/0.5": "0.0892003",
    "relative_suppression_at_node": "0.9108",
    "pair_integral_xi_flat": "0.333332429717",
    "pair_integral_xiprime_flat": "0.141615041677"
  },
  "caveats": [
    "paper decimals are certified lower bounds (truncated); computed values sit slightly above, as expected",
    "R-formula per E-001; D1 per Zeta23/XiPrime/Certificate/D1.lean"
  ],
  "elapsed_seconds": 248.0,
  "timestamp_utc": "2026-08-18T07:50:03.551107+00:00",
  "row_count": 3
}

"""E-013 — Xi-prime certificate constants from the D1 diagonal law (T4).

The paper's xi'-certificate (Remark 7.1) uses the same R-functional as xi
(E-001) but with the pair kernel |u-v| replaced by the xi' diagonal law
(zeta-23-lean Zeta23/XiPrime/Certificate/D1.lean):

    D1(s) = s - 4 s^2 + sum_{k>=0} c_k s^{2k+3},  c_k = 2*4^{k+1} k!/(2k+2)!

    R_xi'(psi) = [∫ψ² + ∬ D1(|u-v|) ψ(u)ψ(v) du dv] / (∫ψ)²   on [-1/2, 1/2]

Mechanism being validated (D-022): the xi' coefficients
C(N;L_T) = -Lambda(N) + (Lambda*log)(N)/L_T + O(l^-2), L_T = l/2 + i pi/4,
Re(1/L_T) ~ 2/l, give the diagonal law a PERFECT-SQUARE node:
  D1(s) = s - 4s^2 + 4s^3 - ... ~= s(1-2s)^2  (leading three terms),
i.e. the -Lambda and +(2/l)(Lambda*log) diagonals CANCEL at prime scale
p ~ e^{l/2} = sqrt(T/2pi). The pair integral for the flat window drops
from 1/3 (xi) to ~0.1416 (xi') — exactly the certificate gain
0.85838 - 0.66667.

Checks:
  V1. flat window:  2 - R_xi'(1) matches certified 0.85838 (within 1e-3;
      paper decimals are certified lower bounds, truncation expected)
  V2. quartic window v(s) = 1 - (7/100)(2s)^2 - (51/200)(2s)^4:
      2 - R_xi' matches certified 0.86864 (within 1e-3)
  V3. mechanism: (a) D1(x)/x -> 1 as x -> 0; (b) deep mid-scale
      suppression D1(0.5)/0.5 < 0.1 (the (2s-1)^2 node, filled to ~0.09
      by higher-order terms); (c) pair-integral reduction equals the
      certificate gain (bookkeeping identity)
  V4. xi sanity: same code path with kernel x reproduces R = 4/3 (flat),
      quad-precision tolerance 1e-5

Outputs: research/numerics/data/E-013-xiprime-R-constants.csv (+ .json)
"""
import json
import time
from datetime import datetime, timezone
from pathlib import Path

import mpmath as mp

REPO = Path(__file__).resolve().parents[3]
DATA = REPO / "research" / "numerics" / "data"
CSV = DATA / "E-013-xiprime-R-constants.csv"
JSON = DATA / "E-013-xiprime-R-constants.json"

K_TAIL = 40          # series terms in D1
CERT_FLAT = 0.85838
CERT_QUARTIC = 0.86864


def d1_coeff(k):
    return 2 * 4 ** (k + 1) * mp.factorial(k) / mp.factorial(2 * k + 2)


def d1(x):
    s = x - 4 * x * x
    xp = x ** 3
    for k in range(K_TAIL):
        s += d1_coeff(k) * xp
        xp *= x * x
    return s


def R(psi, kernel):
    """R(psi) = [∫ψ² + ∬ kernel(|u-v|)ψ(u)ψ(v)] / (∫ψ)² on [-1/2,1/2]."""
    i1 = mp.quad(lambda u: psi(u) ** 2, [-mp.mpf("0.5"), mp.mpf("0.5")])
    i2 = mp.quad(lambda u: mp.quad(
        lambda v: kernel(abs(u - v)) * psi(u) * psi(v),
        [-mp.mpf("0.5"), mp.mpf("0.5")]), [-mp.mpf("0.5"), mp.mpf("0.5")])
    n = mp.quad(psi, [-mp.mpf("0.5"), mp.mpf("0.5")])
    return (i1 + i2) / n ** 2


def flat(u):
    return mp.mpf(1)


def quartic(u):
    x = 2 * u
    return 1 - mp.mpf("0.07") * x ** 2 - mp.mpf("0.255") * x ** 4


def main() -> int:
    import sys
    if CSV.exists() and "--force" not in sys.argv:
        print(f"cached: {CSV} (use --force to recompute)")
        return 0
    t0 = time.time()
    mp.mp.dps = 40

    # V4: xi sanity (kernel x); mp.quad accuracy ~1e-6 here
    r_xi_flat = R(flat, lambda x: x)
    check_v4 = bool(abs(r_xi_flat - mp.mpf(4) / 3) < 1e-5)

    # V1/V2: xi' constants
    r_flat = R(flat, d1)
    r_quartic = R(quartic, d1)
    c_flat = 2 - r_flat
    c_quartic = 2 - r_quartic
    check_v1 = bool(c_flat >= CERT_FLAT and c_flat - CERT_FLAT < 1e-3)
    check_v2 = bool(c_quartic >= CERT_QUARTIC and c_quartic - CERT_QUARTIC < 1e-3)

    # V3: mechanism — (a) D1(x)/x -> 1 as x -> 0; (b) mid-scale node;
    # (c) pair-integral reduction == certificate gain
    ratio_small = d1(mp.mpf("0.001")) / mp.mpf("0.001")
    v3a = bool(abs(ratio_small - (1 - 4 * mp.mpf("0.001"))) < 1e-4)
    node_ratio = d1(mp.mpf("0.5")) / mp.mpf("0.5")
    v3b = bool(node_ratio < 0.1)
    pair_xi = r_xi_flat - 1            # = 1/3
    pair_xip = r_flat - 1
    v3c = bool(abs((pair_xi - pair_xip) - (c_flat - mp.mpf(2) / 3)) < 1e-5)
    check_v3 = bool(v3a and v3b and v3c)
    sup_deficit = 1 - node_ratio       # deepest relative suppression (at node)

    checks = {
        "V1_flat_matches_085838": check_v1,
        "V2_quartic_matches_086864": check_v2,
        "V3_perfect_square_node_mechanism": check_v3,
        "V4_xi_kernel_sanity_4_3": check_v4,
    }
    with CSV.open("w", encoding="utf-8", newline="") as f:
        f.write("window,R_xiprime,certificate_2_minus_R,certified_decimal\n")
        f.write(f"flat,{mp.nstr(r_flat, 15)},{mp.nstr(c_flat, 15)},{CERT_FLAT}\n")
        f.write(f"quartic,{mp.nstr(r_quartic, 15)},{mp.nstr(c_quartic, 15)},"
                f"{CERT_QUARTIC}\n")
        f.write(f"flat_xi_kernel_sanity,{mp.nstr(r_xi_flat, 15)},"
                f"{mp.nstr(2 - r_xi_flat, 15)},0.666666666666667\n")
    meta = {
        "experiment": "E-013",
        "task": "T4: xi'-analogue — why is the method stronger on xi'? (D-022)",
        "description": "R-functional with xi' diagonal law D1 (pair kernel "
                       "x - 4x^2 + ...) reproduces Remark 7.1 constants; "
                       "D1(x) <= x quantifies tight-pair suppression",
        "script": "research/numerics/scripts/E-013-xiprime-R-constants.py",
        "parameters": {"dps": 40, "d1_series_terms": K_TAIL},
        "checks": checks,
        "key_values": {
            "R_xiprime_flat": mp.nstr(r_flat, 15),
            "cert_flat_computed": mp.nstr(c_flat, 15),
            "cert_flat_paper": CERT_FLAT,
            "R_xiprime_quartic": mp.nstr(r_quartic, 15),
            "cert_quartic_computed": mp.nstr(c_quartic, 15),
            "cert_quartic_paper": CERT_QUARTIC,
            "R_xi_flat_sanity": mp.nstr(r_xi_flat, 15),
            "node_ratio_D1(0.5)/0.5": mp.nstr(node_ratio, 6),
            "relative_suppression_at_node": mp.nstr(sup_deficit, 6),
            "pair_integral_xi_flat": mp.nstr(pair_xi, 12),
            "pair_integral_xiprime_flat": mp.nstr(pair_xip, 12),
        },
        "caveats": [
            "paper decimals are certified lower bounds (truncated); "
            "computed values sit slightly above, as expected",
            "R-formula per E-001; D1 per Zeta23/XiPrime/Certificate/D1.lean",
        ],
        "elapsed_seconds": round(time.time() - t0, 1),
        "timestamp_utc": datetime.now(timezone.utc).isoformat(),
        "row_count": 3,
    }
    JSON.write_text(json.dumps(meta, indent=2), encoding="utf-8")
    ok = all(checks.values())
    print(f"xi flat sanity R={mp.nstr(r_xi_flat, 10)} (want 4/3)")
    print(f"xi' flat:    R={mp.nstr(r_flat, 10)}  2-R={mp.nstr(c_flat, 8)} "
          f"(certified {CERT_FLAT})")
    print(f"xi' quartic: R={mp.nstr(r_quartic, 10)}  2-R={mp.nstr(c_quartic, 8)} "
          f"(certified {CERT_QUARTIC})")
    print(f"mechanism: max relative pair suppression on [0,1] = "
          f"{mp.nstr(sup_deficit, 4)}")
    print(f"wrote {CSV.name}, {JSON.name} — "
          f"{'ALL CHECKS PASS' if ok else 'CHECK FAILURE: ' + str(checks)}")
    return 0 if ok else 1


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