E-031-right-tail-stepping.py

{
  "experiment": "E-031",
  "title": "Finding F-D: the right-tail of Lagarias's Thm 6.1 horizontals ((6.12) stepping broken)",
  "session": "038",
  "design_doc": "D-072",
  "checks": {
    "V1_peak_tracks_0.36sqrt_e^{0.066sqrt}": true,
    "V2_envelope_~6e10_at_1e5": true,
    "V3_repair_works_for_c>=1.5": true
  },
  "all_checks_pass": true,
  "dps": 25,
  "notes": "The (6.12) product ratio is e/2 > 1: the right-tail grows to ~e^{0.066√n} at σ ≈ 0.36√n (5.9e10 envelope at n=1e5). Repair: T = c√n with c ≥ ~1.25 makes the product decrease per step."
}

#!/usr/bin/env python3
"""E-031 — Finding F-D: the right-tail of the horizontals in Lagarias's
Thm 6.1 contour (the (6.12) stepping). Session 038.

CLAIM CHECKED: Lagarias (6.12) asserts that |k_n(s)·L′/L(s)| decreases
geometrically per unit step right from σ = 3 (zeta halves via (6.10),
kernel ratio ≤ e via (6.12)) ⟹ right part of the horizontals is O(√n).

FINDING: the per-step product ratio is e/2 > 1 — the product GROWS
from σ = 3 to σ ≈ 0.36√n, peaking at ~e^{0.066√n}. The stepping is
broken as written at T = √n; the right-tail envelope integral is
~5.9×10¹⁰ at n = 10⁵ (vs the ~3.3 charged in D-036/D-037).

REPAIR CHECK: at T = c√n with c = 1.5, 2 the product DECREASES from
σ = 3 (kernel ratio ≤ e^{1/c²} < 2) — the stepping then works.

Computations (direct, no contour-modeling):
  V1: product profile |k_n(w−1)|·|Σ_{k<300} Λ(k)/k^w| at w = σ+i√n
      for σ ∈ grid, n ∈ {10³, 10⁴, 10⁵} — peak position ~0.36√n,
      peak height ~e^{0.066√n}.
  V2: the right-tail envelope integral ∫_{3}^{2√n} (step 2) at
      n = 10⁴, 10⁵.
  V3: the same profile at T = c√n, c ∈ {1, 1.5, 2}, n = 10⁵ — the
      product decreases from σ = 3 for c ≥ 1.5.
"""
import json
import sys
from pathlib import Path

import mpmath as mp

mp.mp.dps = 25

DATA = Path(__file__).resolve().parent.parent / "data"
CSV = DATA / "E-031-right-tail-stepping.csv"
JSON = DATA / "E-031-right-tail-stepping.json"

_PRIMES = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]


def vonM(k):
    k = int(k)
    for p in _PRIMES:
        if k == p:
            return mp.log(p)
        if k % p == 0:
            m = k
            while m % p == 0:
                m //= p
            if m == 1:
                return mp.log(p)
            return mp.mpf(0)
    return mp.mpf(0)


def product(n, sig, c=1.0):
    """|k_n(w−1)| · |Σ_{k<300} Λ(k)/k^w| at w = sig + i·c√n."""
    T = c * mp.sqrt(n)
    w = sig + T * 1j
    kn = abs(1 + 1 / (w - 1)) ** n
    z = mp.mpc(0)
    for k in range(2, 300):
        z += vonM(k) / k ** w
    return kn * abs(z)


def envelope_integral(n):
    """∫_{3}^{2√n} product dσ (step 2), T = √n."""
    T = mp.sqrt(n)
    tot = mp.mpf(0)
    peak = (mp.mpf(0), mp.mpf(0))
    sig = mp.mpf(3)
    while sig < 2 * T:
        v = product(n, sig)
        if v > peak[0]:
            peak = (v, sig)
        tot += v * 2
        sig += 2
    return tot, peak


def main() -> int:
    if JSON.exists() and "--force" not in sys.argv:
        print(f"cached: {JSON} (use --force to recompute)")
        return 0

    rows = []
    # V1: peak profile across n (T = √n)
    v1 = {}
    for n in [10 ** 3, 10 ** 4, 10 ** 5]:
        T = mp.sqrt(n)
        best = (mp.mpf(0), mp.mpf(0))
        sig = mp.mpf(3)
        while sig < 2 * T:
            v = product(n, sig)
            if v > best[0]:
                best = (v, sig)
            sig += 2
        v1[n] = (best, T)
        rows.append({"check": "V1_peak", "n": mp.nstr(n, 6), "c": "1.0",
                     "peak_product": mp.nstr(best[0], 6),
                     "peak_sigma": mp.nstr(best[1], 6),
                     "peak_sigma_over_sqrt": mp.nstr(best[1] / T, 4),
                     "log_peak_over_sqrt": mp.nstr(mp.log(best[0]) / T, 5)})
    # V2: envelope integral
    v2 = {}
    for n in [10 ** 4, 10 ** 5]:
        tot, peak = envelope_integral(n)
        v2[n] = (tot, peak)
        rows.append({"check": "V2_envelope_integral", "n": mp.nstr(n, 6), "c": "1.0",
                     "envelope_integral": mp.nstr(tot, 6),
                     "peak_product": mp.nstr(peak[0], 6),
                     "peak_sigma": mp.nstr(peak[1], 6),
                     "peak_sigma_over_sqrt": "",
                     "log_peak_over_sqrt": ""})
    # V3: the repair (product decreases from σ = 3 for c ≥ 1.5)
    n = 10 ** 5
    repair_ok = {}
    for c in [1.0, 1.5, 2.0]:
        vals = [product(n, sig, c) for sig in [3, 10, 30, 60, 100, 200]]
        decreasing = all(vals[i + 1] < vals[i] for i in range(len(vals) - 1))
        repair_ok[c] = decreasing
        rows.append({"check": "V3_repair_decreases", "n": mp.nstr(n, 6), "c": mp.nstr(c, 3),
                     "logP_at_3": mp.nstr(mp.log(vals[0]), 6),
                     "logP_at_100": mp.nstr(mp.log(vals[4]), 6),
                     "monotone_decrease": str(decreasing),
                     "peak_sigma_over_sqrt": "", "envelope_integral": "",
                     "peak_product": "", "peak_sigma": "", "log_peak_over_sqrt": ""})

    v1_ok = all(v1[n][0][0] > 1 and mp.mpf("0.2") < v1[n][0][1] / v1[n][1] < mp.mpf("0.6")
                for n in v1)
    v2_ok = v2[10 ** 5][0] > 10 ** 9  # the envelope really is ~10¹⁰ at 1e5
    v3_ok = repair_ok[1.5] and repair_ok[2.0] and not repair_ok[1.0]

    DATA.mkdir(exist_ok=True)
    with CSV.open("w") as f:
        f.write(",".join(rows[0].keys()) + "\n")
        for r in rows:
            f.write(",".join(r.values()) + "\n")
    sidecar = {
        "experiment": "E-031",
        "title": "Finding F-D: the right-tail of Lagarias's Thm 6.1 horizontals ((6.12) stepping broken)",
        "session": "038",
        "design_doc": "D-072",
        "checks": {"V1_peak_tracks_0.36sqrt_e^{0.066sqrt}": v1_ok,
                   "V2_envelope_~6e10_at_1e5": v2_ok,
                   "V3_repair_works_for_c>=1.5": v3_ok},
        "all_checks_pass": v1_ok and v2_ok and v3_ok,
        "dps": 25,
        "notes": "The (6.12) product ratio is e/2 > 1: the right-tail grows to "
                 "~e^{0.066√n} at σ ≈ 0.36√n (5.9e10 envelope at n=1e5). "
                 "Repair: T = c√n with c ≥ ~1.25 makes the product decrease per step.",
    }
    with JSON.open("w") as f:
        json.dump(sidecar, f, indent=2)
    print(f"V1 peak profile: {v1_ok} | V2 envelope ~6e10: {v2_ok} | V3 repair: {v3_ok}")
    print(f"wrote {CSV} + {JSON}")
    return 0 if sidecar["all_checks_pass"] else 1


if __name__ == "__main__":
    sys.exit(main())