#!/usr/bin/env python3
"""
PROVENANCE: NUMERICAL

Computation 202 -- the Higgs self-coupling test at three loops on current inputs:
three-loop running is worth nothing at M_*, the 2013 anchor is not reproducible,
and the honest agreement is 1.6 sigma on the 2024 Review and 1.4 sigma on the 2026 one
==============================================================================
WHAT THIS DOES.  The book's stated test of the bridge premise is comp 123's
lambda_SM(M_*) = 0.0927 +- 0.0007 against the substrate's e^-1/4 = 0.09197, a 1.0
sigma agreement anchored to a literature value rather than computed.  This file
computes the number itself, at full three loops, with the NNLO matching conditions
at the top mass, on the inputs the anchor used and on the current ones.

SOURCES.  Matching at mu = M_t: Buttazzo et al. 2013 (arXiv:1307.3536) eqs. 55 and
57-60, which are their NNLO fits to M_h, M_t, alpha_3(M_Z) and M_W.  Running: the
same paper's appendix B, the complete three-loop SM beta functions for g1 (GUT
normalised), g2, g3 (with the four-loop pure QCD term), y_t and lambda, with y_b and
y_tau at two and one loops.  Inputs: Buttazzo 2013 (M_t 173.34 +- 0.76, M_h
125.15 +- 0.24, alpha_3 0.1184 +- 0.0007, M_W 80.384 +- 0.014) and PDG 2024 (M_t
172.57 +- 0.29 direct, M_h 125.20 +- 0.11, alpha_3 0.1180 +- 0.0009, M_W
80.3692 +- 0.0133), plus the pole-mass ambiguity of order 0.5 GeV the PDG top review
states between the direct-measurement mass and the pole mass the matching needs.
Theory errors: +- 0.00030 on lambda(M_t) and +- 0.00050 on y_t(M_t), the paper's.

RESULTS.
  1. VALIDATION.  Run to the Planck mass on the 2013 central inputs, the five
     couplings reproduce the paper's Planck-scale values (g1 0.6154, g2 0.5055,
     g3 0.4873, y_t 0.3825, lambda -0.0143) to four decimals, so the three-loop
     system is coded right.  comp 198's two-loop 0.09228 is reproduced on its inputs.
  2. THREE LOOPS ARE WORTH NOTHING HERE.  The third loop moves lambda(M_*) by
     +0.00007, a tenth of the error.  The 0.0004 that comp 198 attributed to three-loop
     running and threshold matching does not exist on the running side.
  3. THE ANCHOR IS NOT REPRODUCIBLE.  On Buttazzo's own inputs their own formulas give
     0.0924 +- 0.0011 at every loop order, not 0.0927.  comp 123's 1.0 sigma was
     0.4 sigma on 2013 inputs, with an error budget (0.0007) that understated the
     2013 top-mass error.
  4. CURRENT INPUTS.  The top mass fell from 173.34 to 172.57 GeV and a lighter top
     raises lambda at M_*: PDG 2024 inputs give lambda(M_*) = 0.0932 +- 0.0006 from
     the quoted errors alone (2.2 sigma) and 0.0932 +- 0.0008 once the pole-mass
     ambiguity is counted (1.6 sigma); the cross-section pole mass 172.4 +- 0.7 gives
     0.0934 +- 0.0009, also 1.6 sigma.  The honest statement on those inputs is 1.6
     sigma, 1.4 percent, and it is what the book printed from v32.137 to v32.152.
  5. WHERE EXACTNESS WOULD SIT.  lambda(M_*) = e^-1/4 exactly needs M_t = 173.7 GeV,
     four sigma above the direct measurement and two above the cross-section mass.
  6. THE SAME NUMBER IN c (comp 176).  Re-solving the bridge premise for the kinetic
     coefficient on 2024 inputs moves c off 1 by the same 1.6 sigma, and sqrt(e) stays
     excluded by more than a hundred sigma; nothing about the reading-B exclusion
     changes.
  8. THE 2026 REVIEW (v32.153 re-check).  PDG 2026 moves the inputs to m_t 172.60 +-
     0.27, m_h 125.13 +- 0.11, M_W 80.3625 +- 0.0077 with alpha_s unchanged at 0.1180 +-
     0.0009.  The Higgs mass's 0.07 GeV does most of the work: lambda(M_*) = 0.0930 +-
     0.0008 with the pole-mass ambiguity counted, 1.4 sigma and 1.2 percent above
     e^-1/4 (1.9 sigma on the quoted errors alone).  The book prints 1.4 sigma from
     v32.153, on the cover too; the recurring task is to re-run this part at each
     edition of the Review.

TAMPERS.  Dropping the three-loop QCD piece of the y_t matching (0.94018 instead of
0.93690) moves lambda(M_*) by -0.0006, nine times the three-loop running step, so the
precision lives in the matching, not in the running; and the input set that reproduces
comp 123's 0.0007 budget (m_t +- 0.5, m_h +- 0.17) is shown to be neither 2013's nor
2024's.

SCOPE.  This file supersedes comp 123's number and the loop-order sequence comp 198
quoted; both files carry scope notes pointing here.  It does not touch the substrate
side: e^-1/4 is comp 100's and comp 174's, exact and D-independent.

Runs in a few seconds with the standard library only.
"""

import math

PI2 = 16 * math.pi ** 2
MPL = 1.2209e19
PST = 0.25 * math.exp(-1.0)

# ---------------- matching at mu = M_t (Buttazzo et al. 2013) ----------------
def matching(Mt=173.34, Mh=125.15, a3=0.1184, MW=80.384, dlam_th=0.0, dyt_th=0.0, yt_nnlo_only=False):
    das = (a3 - 0.1184) / 0.0007
    dmw = (MW - 80.384) / 0.014
    lam = 0.12604 + 0.00206 * (Mh - 125.15) - 0.00004 * (Mt - 173.34) + dlam_th
    yt = (0.94018 if yt_nnlo_only else 0.93690) + 0.00556 * (Mt - 173.34) - 0.00042 * das + dyt_th
    g2 = 0.64779 + 0.00004 * (Mt - 173.34) + 0.00011 * dmw
    gY = 0.35830 + 0.00011 * (Mt - 173.34) - 0.00020 * dmw
    g3 = 1.1666 + 0.00314 * das - 0.00046 * (Mt - 173.34)
    g1 = math.sqrt(5.0 / 3.0) * gY
    yb = math.sqrt(2) * 2.75 / 246.22      # m_b(M_t) ~ 2.75 GeV from m_b(m_b) = 4.2 GeV
    ytau = math.sqrt(2) * 1.777 / 246.22
    return dict(g1=g1, g2=g2, g3=g3, yt=yt, lam=lam, yb=yb, ytau=ytau)

# ---------------- RGEs, d/d ln mu^2 of squared couplings (appendix B) ----------------
def rhs(Y, loops=3, light=True):
    G1, G2, G3, T, lam, B, Tau = Y          # squares: g1^2, g2^2, g3^2, yt^2, lambda, yb^2, ytau^2
    if not light:
        B = Tau = 0.0
    k1, k2, k3 = 1 / (4 * math.pi) ** 2, 1 / (4 * math.pi) ** 4, 1 / (4 * math.pi) ** 6
    # gauge
    d1 = G1 ** 2 * k1 * (41 / 10)
    d2 = G2 ** 2 * k1 * (-19 / 6)
    d3 = G3 ** 2 * k1 * (-7)
    if loops >= 2:
        d1 += G1 ** 2 * k2 * (44 * G3 / 5 + 27 * G2 / 10 + 199 * G1 / 50 - 17 * T / 10 - B / 2 - 3 * Tau / 2)
        d2 += G2 ** 2 * k2 * (12 * G3 + 35 * G2 / 6 + 9 * G1 / 10 - 3 * T / 2 - 3 * B / 2 - Tau / 2)
        d3 += G3 ** 2 * k2 * (-26 * G3 + 9 * G2 / 2 + 11 * G1 / 10 - 2 * T - 2 * B)
    if loops >= 3:
        d1 += G1 ** 2 * k3 * (T * (189 * T / 16 - 29 * G3 / 5 - 471 * G2 / 32 - 2827 * G1 / 800)
                              + lam * (-9 * lam / 5 + 9 * G2 / 10 + 27 * G1 / 50)
                              + 297 * G3 ** 2 / 5 + 789 * G2 ** 2 / 64 - 388613 * G1 ** 2 / 24000
                              - 3 * G3 * G2 / 5 - 137 * G3 * G1 / 75 + 123 * G2 * G1 / 160)
        d2 += G2 ** 2 * k3 * (T * (147 * T / 16 - 7 * G3 - 729 * G2 / 32 - 593 * G1 / 160)
                              + lam * (-3 * lam + 3 * G2 / 2 + 3 * G1 / 10)
                              + 81 * G3 ** 2 + 324953 * G2 ** 2 / 1728 - 5597 * G1 ** 2 / 1600
                              + 39 * G3 * G2 - G3 * G1 / 5 + 873 * G2 * G1 / 160)
        d3 += G3 ** 2 * k3 * (T * (15 * T - 40 * G3 - 93 * G2 / 8 - 101 * G1 / 40)
                              + 65 * G3 ** 2 / 2 + 109 * G2 ** 2 / 8 - 523 * G1 ** 2 / 120
                              + 21 * G3 * G2 + 77 * G3 * G1 / 15 - 3 * G2 * G1 / 40)
        d3 += G3 ** 5 / (4 * math.pi) ** 8 * (-2472.28)
    # top Yukawa (squared)
    dT = T * k1 * (9 * T / 2 + 3 * B / 2 + Tau - 8 * G3 - 9 * G2 / 4 - 17 * G1 / 20)
    if loops >= 2:
        dT += T * k2 * (T * (-12 * T - 11 * B / 4 - 9 * Tau / 4 - 12 * lam + 36 * G3 + 225 * G2 / 16 + 393 * G1 / 80)
                        + B * (-B / 4 + 5 * Tau / 4 + 4 * G3 + 99 * G2 / 16 + 7 * G1 / 80)
                        + Tau * (-9 * Tau / 4 + 15 * G2 / 8 + 15 * G1 / 8)
                        + 6 * lam ** 2 - 108 * G3 ** 2 - 23 * G2 ** 2 / 4 + 1187 * G1 ** 2 / 600
                        + 9 * G3 * G2 + 19 * G3 * G1 / 15 - 9 * G2 * G1 / 20)
    if loops >= 3:
        dT += T * k3 * (T ** 2 * (58.6028 * T + 198 * lam - 157 * G3 - 1593 * G2 / 16 - 2437 * G1 / 80)
                        + lam * T * (15 * lam / 4 + 16 * G3 - 135 * G2 / 2 - 127 * G1 / 10)
                        + T * (363.764 * G3 ** 2 + 16.990 * G2 ** 2 - 24.422 * G1 ** 2 + 48.370 * G3 * G2
                               + 18.074 * G3 * G1 + 34.829 * G2 * G1)
                        + lam ** 2 * (-36 * lam + 45 * G2 + 9 * G1)
                        + lam * (-171 * G2 ** 2 / 16 - 1089 * G1 ** 2 / 400 + 117 * G2 * G1 / 40)
                        - 619.35 * G3 ** 3 + 169.829 * G2 ** 3 + 16.099 * G1 ** 3 + 73.654 * G3 ** 2 * G2
                        - 15.096 * G3 ** 2 * G1 - 21.072 * G3 * G2 ** 2 - 22.319 * G3 * G1 ** 2
                        - 321 / 20 * G3 * G2 * G1 - 4.743 * G2 ** 2 * G1 - 4.442 * G2 * G1 ** 2)
    # quartic
    dl = k1 * (lam * (12 * lam + 6 * T + 6 * B + 2 * Tau - 9 * G2 / 2 - 9 * G1 / 10)
               - 3 * T ** 2 - 3 * B ** 2 - Tau ** 2 + 9 * G2 ** 2 / 16 + 27 * G1 ** 2 / 400 + 9 * G2 * G1 / 40)
    if loops >= 2:
        dl += k2 * (lam ** 2 * (-156 * lam - 72 * T - 72 * B - 24 * Tau + 54 * G2 + 54 * G1 / 5)
                    + lam * T * (-3 * T / 2 - 21 * B + 40 * G3 + 45 * G2 / 4 + 17 * G1 / 4)
                    + lam * B * (-3 * B / 2 + 40 * G3 + 45 * G2 / 4 + 5 * G1 / 4)
                    + lam * Tau * (-Tau / 2 + 15 * G2 / 4 + 15 * G1 / 4)
                    + lam * (-73 * G2 ** 2 / 16 + 1887 * G1 ** 2 / 400 + 117 * G2 * G1 / 40)
                    + T ** 2 * (15 * T - 3 * B - 16 * G3 - 4 * G1 / 5)
                    + T * (-9 * G2 ** 2 / 8 - 171 * G1 ** 2 / 200 + 63 * G2 * G1 / 20)
                    + B ** 2 * (-3 * T + 15 * B - 16 * G3 + 2 * G1 / 5)
                    + B * (-9 * G2 ** 2 / 8 + 9 * G1 ** 2 / 40 + 27 * G2 * G1 / 20)
                    + Tau ** 2 * (5 * Tau - 6 * G1 / 5) + Tau * (-3 * G2 ** 2 / 8 - 9 * G1 ** 2 / 8 + 33 * G2 * G1 / 20)
                    + 305 * G2 ** 3 / 32 - 3411 * G1 ** 3 / 4000 - 289 * G2 ** 2 * G1 / 160 - 1677 * G2 * G1 ** 2 / 800)
    if loops >= 3:
        dl += k3 * (lam ** 3 * (6011.35 * lam + 873 * T - 387.452 * G2 - 77.490 * G1)
                    + lam ** 2 * T * (1768.26 * T + 160.77 * G3 - 359.539 * G2 - 63.869 * G1)
                    + lam ** 2 * (-790.28 * G2 ** 2 - 185.532 * G1 ** 2 - 316.64 * G2 * G1)
                    + lam * T ** 2 * (-223.382 * T - 662.866 * G3 - 5.470 * G2 - 21.015 * G1)
                    + lam * T * (356.968 * G3 ** 2 - 319.664 * G2 ** 2 - 74.8599 * G1 ** 2 + 15.1443 * G3 * G2
                                 + 17.454 * G3 * G1 + 5.615 * G2 * G1)
                    + lam * G2 ** 2 * (-57.144 * G3 + 865.483 * G2 + 79.638 * G1)
                    + lam * G1 ** 2 * (-8.381 * G3 + 61.753 * G2 + 28.168 * G1)
                    + T ** 3 * (-243.149 * T + 250.494 * G3 + 74.138 * G2 + 33.930 * G1)
                    + T ** 2 * (-50.201 * G3 ** 2 + 15.884 * G2 ** 2 + 15.948 * G1 ** 2 + 13.349 * G3 * G2
                                + 17.570 * G3 * G1 - 70.356 * G2 * G1)
                    + T * G3 * (16.464 * G2 ** 2 + 1.016 * G1 ** 2 + 11.386 * G2 * G1)
                    + T * G2 ** 2 * (62.500 * G2 + 13.041 * G1)
                    + T * G1 ** 2 * (10.627 * G2 + 11.117 * G1)
                    + G3 * (7.536 * G2 ** 3 + 0.663 * G1 ** 3 + 1.507 * G2 ** 2 * G1 + 1.105 * G2 * G1 ** 2)
                    - 114.091 * G2 ** 4 - 1.508 * G1 ** 4 - 37.889 * G2 ** 3 * G1 + 6.500 * G2 ** 2 * G1 ** 2
                    - 1.543 * G2 * G1 ** 3)
    # bottom and tau Yukawa (squared), two and one loop
    dB = B * k1 * (3 * T / 2 + 9 * B / 2 + Tau - 8 * G3 - 9 * G2 / 4 - G1 / 4)
    if loops >= 2:
        dB += B * k2 * (T * (-T / 4 - 11 * B / 4 + 5 * Tau / 4 + 4 * G3 + 99 * G2 / 16 + 91 * G1 / 80)
                        + B * (-12 * B - 9 * Tau / 4 - 12 * lam + 36 * G3 + 225 * G2 / 16 + 237 * G1 / 80)
                        + Tau * (-9 * Tau / 4 + 15 * G2 / 8 + 15 * G1 / 8)
                        + 6 * lam ** 2 - 108 * G3 ** 2 - 23 * G2 ** 2 / 4 - 127 * G1 ** 2 / 600
                        + 9 * G3 * G2 + 31 * G3 * G1 / 15 - 27 * G2 * G1 / 20)
    dTau = Tau * k1 * (5 * Tau / 2 + 3 * T + 3 * B - 9 * G2 / 4 - 27 * G1 / 20)
    if not light:
        dB = dTau = 0.0
    return [d1, d2, d3, dT, dl, dB, dTau]

def run(bc, mu0, mu1, loops=3, light=True, n=2000):
    Y = [bc["g1"] ** 2, bc["g2"] ** 2, bc["g3"] ** 2, bc["yt"] ** 2, bc["lam"], bc["yb"] ** 2, bc["ytau"] ** 2]
    t0, t1 = 2 * math.log(mu0), 2 * math.log(mu1)
    h = (t1 - t0) / n
    for _ in range(n):
        k1 = rhs(Y, loops, light)
        k2 = rhs([y + 0.5 * h * k for y, k in zip(Y, k1)], loops, light)
        k3 = rhs([y + 0.5 * h * k for y, k in zip(Y, k2)], loops, light)
        k4 = rhs([y + h * k for y, k in zip(Y, k3)], loops, light)
        Y = [y + h / 6 * (a + 2 * b + 2 * c + d) for y, a, b, c, d in zip(Y, k1, k2, k3, k4)]
    return dict(g1=math.sqrt(Y[0]), g2=math.sqrt(Y[1]), g3=math.sqrt(Y[2]), yt=math.sqrt(Y[3]), lam=Y[4])

def mstar(Mh):
    return 4 * math.pi * Mh * math.sqrt(2.0 / 3.0)

def lam_at_mstar(loops=3, light=True, **kw):
    bc = matching(**kw)
    Mh = kw.get("Mh", 125.15)
    return run(bc, kw.get("Mt", 173.34), mstar(Mh), loops, light, n=400)["lam"]

IN13 = dict(Mt=173.34, Mh=125.15, a3=0.1184, MW=80.384)
ER13 = dict(Mt=0.76, Mh=0.24, a3=0.0007, MW=0.014)
IN24 = dict(Mt=172.57, Mh=125.20, a3=0.1180, MW=80.3692)
ER24 = dict(Mt=0.29, Mh=0.11, a3=0.0009, MW=0.0133)
IN26 = dict(Mt=172.60, Mh=125.13, a3=0.1180, MW=80.3625)   # PDG 2026 (v32.153 re-check)
ER26 = dict(Mt=0.27, Mh=0.11, a3=0.0009, MW=0.0077)
POLE_AMBIGUITY = 0.5           # GeV, PDG top-quark review: direct mass vs pole mass
E = math.e


def budget(loops, inputs, errs, name, quiet=False):
    """Central value and quadrature error of lambda(M_*) from one-sigma shifts of each
    input plus the two matching theory errors."""
    c = lam_at_mstar(loops, **inputs)
    parts = {}
    for k, e in errs.items():
        up = dict(inputs); up[k] = inputs[k] + e
        parts[k] = lam_at_mstar(loops, **up) - c
    parts["lam_th"] = lam_at_mstar(loops, dlam_th=0.00030, **inputs) - c
    parts["yt_th"] = lam_at_mstar(loops, dyt_th=0.00050, **inputs) - c
    tot = math.sqrt(sum(v * v for v in parts.values()))
    sig = (c - PST) / tot
    if not quiet:
        print(f"      {name}: lambda(M_*) = {c:.5f} +- {tot:.5f}  ->  {sig:+.2f} sigma from e^-1/4, {(c - PST) / PST * 100:+.2f} %")
        print("          budget: " + ", ".join(f"{k} {v:+.5f}" for k, v in parts.items()))
    return c, tot, sig


def solve_c(loops, inputs, dlam_th=0.0, dyt_th=0.0):
    """comp 176's equation, lambda_SM(c M_*) = e^-1/(4 c^2), on this file's running;
    the residual is monotone on [0.5, 3], so bisection finds the unique root."""
    bc = matching(dlam_th=dlam_th, dyt_th=dyt_th, **inputs)
    base = mstar(inputs["Mh"])
    def res(c):
        return run(bc, inputs["Mt"], base * c, loops, True, n=400)["lam"] - 1 / (4 * E * c * c)
    # The residual is monotone on [0.5, 3]; the Illinois variant of regula falsi reaches the
    # root to 1e-12 in about ten evaluations where the bisection this used until v32.162
    # took a hundred (the two agree to 1e-12, far inside every printed digit).
    lo, hi = 0.5, 3.0
    flo, fhi = res(lo), res(hi)
    if flo * fhi > 0:
        raise ValueError("solve_c: the root is not bracketed on [0.5, 3]")
    side = 0
    for _ in range(100):
        mid = hi - fhi * (hi - lo) / (fhi - flo)
        fmid = res(mid)
        if abs(fmid) < 1e-15 or abs(hi - lo) < 1e-12:
            return mid
        if fmid * fhi < 0:
            lo, flo = hi, fhi; side = 0
        else:
            if side == 1:
                flo *= 0.5
            side = 1
        hi, fhi = mid, fmid
        if abs(hi - lo) < 1e-12:
            return 0.5 * (lo + hi)
    return 0.5 * (lo + hi)


def c_budget(loops, inputs, errs):
    c0 = solve_c(loops, inputs)
    parts = []
    for k, e in errs.items():
        up = dict(inputs); up[k] = inputs[k] + e
        parts.append(solve_c(loops, up) - c0)
    parts.append(solve_c(loops, inputs, dlam_th=0.00030) - c0)
    parts.append(solve_c(loops, inputs, dyt_th=0.00050) - c0)
    return c0, math.sqrt(sum(p * p for p in parts))


def part1():
    print("  PART 1.  Validation")
    ok = True
    r = run(matching(), 173.34, MPL, 3, True, n=4000)
    ref = dict(g1=0.6154, g2=0.5055, g3=0.4873, yt=0.3825, lam=-0.0143)
    print("      three loops to M_Pl on the 2013 central inputs, against Buttazzo et al.'s Planck-scale fits:")
    for k in ref:
        print(f"          {k}: {r[k]:+.5f}  (paper {ref[k]:+.4f})")
        ok &= abs(r[k] - ref[k]) < 1.5e-4
    bc198 = dict(g1=0.46270, g2=0.64779, g3=1.16655, yt=0.93690, lam=0.12604, yb=0.0, ytau=0.0)
    l2 = run(bc198, 173.34, mstar(125.25), 2, False, n=400)["lam"]
    print(f"      comp 198's two-loop run on its own inputs: {l2:.5f} (comp 198: 0.09228)")
    ok &= abs(l2 - 0.09228) < 3e-5
    return ok


def part2():
    print("  PART 2.  What the third loop is worth at M_*, and the loop-order sequence")
    ok = True
    for nm, inp in (("2013", IN13), ("2024", IN24)):
        seq = [lam_at_mstar(L, **inp) for L in (1, 2, 3)]
        print(f"      {nm} inputs: one loop {seq[0]:.5f}, two loops {seq[1]:.5f}, three loops {seq[2]:.5f};  third-loop step {seq[2] - seq[1]:+.6f}")
        ok &= abs(seq[2] - seq[1]) < 1e-4 and seq[1] > seq[0]
    print("      the third loop moves lambda(M_*) by less than 0.0001; the running is converged at two loops for this purpose")
    return ok


def part3():
    print("  PART 3.  The 2013 anchor: Buttazzo's formulas on Buttazzo's inputs")
    c, tot, sig = budget(3, IN13, ER13, "2013 inputs, three loops")
    budget(2, IN13, ER13, "2013 inputs, two loops")
    print(f"      comp 123 anchored to 0.0927 +- 0.0007; the formulas give {c:.4f}: the anchor is {(0.0927 - c) / tot:+.1f} sigma off and is not reproduced at any loop order")
    print(f"      comp 123's 1.0 sigma was {sig:.1f} sigma on these inputs, with an error that counted m_t at +-0.5 GeV where 2013 had +-0.76")
    return abs(c - 0.0927) > 0.0002 and abs(sig) < 0.6


def part4():
    print("  PART 4.  Current inputs (PDG 2024)")
    c_q, t_q, s_q = budget(3, IN24, ER24, "2024 inputs, quoted errors only")
    er = dict(ER24); er["Mt"] = math.sqrt(ER24["Mt"] ** 2 + POLE_AMBIGUITY ** 2)
    c_a, t_a, s_a = budget(3, IN24, er, f"2024 inputs, pole-mass ambiguity {POLE_AMBIGUITY} GeV added to m_t")
    in_x = dict(IN24); in_x["Mt"] = 172.4
    er_x = dict(ER24); er_x["Mt"] = 0.7
    c_x, t_x, s_x = budget(3, in_x, er_x, "2024 cross-section pole mass 172.4 +- 0.7")
    print(f"      the honest statement: lambda_SM(M_*) = {c_a:.4f} +- {t_a:.4f}, {s_a:.1f} sigma and {(c_a - PST) / PST * 100:.1f} percent above e^-1/4, on either definition of the top mass")
    return 2.0 < s_q < 2.5 and 1.5 < s_a < 1.75 and 1.5 < s_x < 1.75


def part5():
    print("  PART 5.  The top mass at which the match would be exact (other 2024 inputs central)")
    lo, hi = 165.0, 180.0
    for _ in range(40):
        mid = 0.5 * (lo + hi)
        kw = dict(IN24); kw["Mt"] = mid
        if lam_at_mstar(3, **kw) > PST: lo = mid
        else: hi = mid
    print(f"      M_t = {lo:.2f} GeV: {(lo - 172.57) / 0.29:+.1f} sigma from the direct 172.57 +- 0.29, {(lo - 172.4) / 0.7:+.1f} sigma from the cross-section 172.4 +- 0.7")
    return 173.5 < lo < 174.0


def part6():
    print("  PART 6.  The same agreement in comp 176's variable, the kinetic coefficient c")
    er = dict(ER24); er["Mt"] = math.sqrt(ER24["Mt"] ** 2 + POLE_AMBIGUITY ** 2)
    c0, dc = c_budget(3, IN24, er)
    c13, dc13 = c_budget(3, IN13, ER13)
    print(f"      2013 inputs: c = {c13:.4f} +- {dc13:.4f}, {(c13 - 1) / dc13:+.1f} sigma from c = 1 (comp 176 on its anchor: 0.9978 +- 0.0054, 0.4 sigma)")
    print(f"      2024 inputs: c = {c0:.4f} +- {dc:.4f}, {(c0 - 1) / dc:+.1f} sigma from c = 1, {(math.sqrt(E) - c0) / dc:.0f} sigma from reading B's sqrt(e) = {math.sqrt(E):.4f}")
    print("      the number is the Higgs test in another variable, one coincidence and not two; reading B stays excluded by more than a hundred sigma")
    return abs(abs((c0 - 1) / dc) - 1.6) < 0.25 and (math.sqrt(E) - c0) / dc > 100


def part7():
    print("  PART 7.  Tampers")
    ok = True
    base = lam_at_mstar(3, **IN13)
    alt = lam_at_mstar(3, yt_nnlo_only=True, **IN13)
    print(f"      y_t matched without its three-loop QCD piece (0.94018 for 0.93690): lambda(M_*) moves by {alt - base:+.5f}, {abs(alt - base) / 0.00007:.0f} times the third-loop running step; the precision lives in the matching")
    ok &= abs(alt - base) > 4e-4
    c, tot, _ = budget(3, IN13, dict(Mt=0.5, Mh=0.17, a3=0.0009, MW=0.014), "comp 123's error set (m_t 0.5, m_h 0.17) on 2013 central values", quiet=True)
    print(f"      comp 123's error set on 2013 central values: {c:.5f} +- {tot:.5f}; the +-0.0007 is reproduced only by shrinking the 2013 top-mass error, and the central value still is not 0.0927")
    ok &= abs(tot - 0.0007) < 0.0002 and abs(c - 0.0927) > 0.0002
    return ok


def part8():
    print("  PART 8.  Re-check on PDG 2026 inputs (m_t 172.60 +- 0.27, m_h 125.13 +- 0.11, alpha_s 0.1180 +- 0.0009, M_W 80.3625 +- 0.0077)")
    c_q, t_q, s_q = budget(3, IN26, ER26, "2026 inputs, quoted errors only")
    er = dict(ER26); er["Mt"] = math.sqrt(ER26["Mt"] ** 2 + POLE_AMBIGUITY ** 2)
    c_a, t_a, s_a = budget(3, IN26, er, f"2026 inputs, pole-mass ambiguity {POLE_AMBIGUITY} GeV added to m_t")
    c24, t24, s24 = budget(3, IN24, dict(ER24, Mt=math.sqrt(ER24["Mt"] ** 2 + POLE_AMBIGUITY ** 2)), "2024 inputs, same treatment", quiet=True)
    print(f"      the move from 2024 to 2026 inputs: lambda_SM(M_*) {c24:.5f} -> {c_a:.5f} ({c_a - c24:+.5f}, the Higgs mass's 0.07 GeV doing most of it), sigma {s24:.2f} -> {s_a:.2f}")
    print(f"      the stated test on current inputs: lambda_SM(M_*) = {c_a:.4f} +- {t_a:.4f}, {s_a:.1f} sigma and {(c_a - PST) / PST * 100:.1f} percent above e^-1/4")
    seq = [lam_at_mstar(L, **IN26) for L in (1, 2, 3)]
    print(f"      loop-order sequence on 2026 inputs: {seq[0]:.5f} -> {seq[1]:.5f} -> {seq[2]:.5f}; Z^2 = lambda/b = {4 * c_a:.3f}, Z^2/e^-1 = {4 * c_a * E:.4f}")
    lo, hi = 165.0, 180.0
    for _ in range(40):
        mid = 0.5 * (lo + hi)
        kw = dict(IN26); kw["Mt"] = mid
        if lam_at_mstar(3, **kw) > PST: lo = mid
        else: hi = mid
    print(f"      exact agreement on 2026 inputs would need M_t = {lo:.2f} GeV, {(lo - 172.60) / 0.27:+.1f} sigma from the direct 172.60 +- 0.27")
    er26 = dict(ER26); er26["Mt"] = math.sqrt(ER26["Mt"] ** 2 + POLE_AMBIGUITY ** 2)
    c0, dc = c_budget(3, IN26, er26)
    print(f"      in comp 176's variable: c = {c0:.4f} +- {dc:.4f}, {(c0 - 1) / dc:+.1f} sigma from c = 1, {(math.sqrt(E) - c0) / dc:.0f} sigma from sqrt(e)")
    return abs(c_a - c24) < 0.0003 and abs(s_a - s24) < 0.3 and 173.0 < lo < 174.0


def run_vec(Y, mu0, mu1, loops=3, light=True, n=2000):
    """run() on the full seven-coupling vector, returning the vector (yb, ytau included) so a
    trajectory can be continued or re-entered with one coupling replaced."""
    t0, t1 = 2 * math.log(mu0), 2 * math.log(mu1)
    h = (t1 - t0) / n
    for _ in range(n):
        k1 = rhs(Y, loops, light)
        k2 = rhs([y + 0.5 * h * k for y, k in zip(Y, k1)], loops, light)
        k3 = rhs([y + 0.5 * h * k for y, k in zip(Y, k2)], loops, light)
        k4 = rhs([y + h * k for y, k in zip(Y, k3)], loops, light)
        Y = [y + h / 6 * (a + 2 * b + 2 * c + d) for y, a, b, c, d in zip(Y, k1, k2, k3, k4)]
    return Y


def part9():
    """The prediction run back down.  The book had said that running e^-1/4 from M_* down to
    the electroweak scale returns lambda(v) = 0.1294 and hence v = m_h / sqrt(2 lambda) =
    246.2 GeV, 0.01 percent from the observed value.  That applied the tree-level relation
    to the MS-bar running coupling at mu = v, which is not the coupling the relation is
    about: on the measured trajectory itself lambda(v) = 0.119 and the same arithmetic gives
    256 GeV.  The honest statement runs the prediction down to the top mass, where the
    NNLO matching ties lambda to (m_h, v), and reads v off the matched relation."""
    print("  PART 9.  The prediction run back down: what e^-1/4 at M_* says about v")
    bc = matching(**IN26)
    Y0 = [bc["g1"] ** 2, bc["g2"] ** 2, bc["g3"] ** 2, bc["yt"] ** 2, bc["lam"], bc["yb"] ** 2, bc["ytau"] ** 2]
    Mt, Mh, v = IN26["Mt"], IN26["Mh"], 246.22
    Yup = run_vec(list(Y0), Mt, mstar(Mh))
    Yv_meas = run_vec(list(Y0), Mt, v)
    Yp = list(Yup); Yp[4] = PST
    lam_mt_pred = run_vec(list(Yp), mstar(Mh), Mt)[4]
    lam_v_pred = run_vec(list(Yp), mstar(Mh), v)[4]
    v_pred = v * math.sqrt(bc["lam"] / lam_mt_pred)      # the matched lambda(m_t) scales as m_h^2 / v^2 at leading order
    print(f"      measured trajectory: lambda(m_t) = {bc['lam']:.5f} (NNLO matching from m_h, v), lambda(M_*) = {Yup[4]:.5f}, lambda(v) = {Yv_meas[4]:.5f}")
    print(f"      the old arithmetic on the measured trajectory: m_h / sqrt(2 lambda(v)) = {Mh / math.sqrt(2 * Yv_meas[4]):.1f} GeV, so it never returned 246.2 even for the data")
    print(f"      e^-1/4 at M_* run down: lambda(m_t) = {lam_mt_pred:.5f} ({(lam_mt_pred / bc['lam'] - 1) * 100:+.2f} percent against the matched value), lambda(v) = {lam_v_pred:.5f}")
    print(f"      v from m_h alone, through the matched relation: {v_pred:.1f} GeV, {(v_pred / v - 1) * 100:+.2f} percent against 246.22")
    print(f"      the book's earlier 0.01 percent is withdrawn; the deficit in v is the 1.2 percent deficit in lambda seen from the other end, halved by the square root")
    return 0.124 < lam_mt_pred < 0.126 and 246.5 < v_pred < 248.0 and Yv_meas[4] < 0.121


def main():
    print("=" * 100)
    print("Computation 202: the Higgs self-coupling test at three loops on current inputs")
    print("=" * 100)
    results = [
        ("three-loop system validated against the paper's Planck-scale values; comp 198 reproduced", part1()),
        ("the third loop moves lambda(M_*) by less than 0.0001", part2()),
        ("the 0.0927 anchor is not reproduced; the 2013 agreement is 0.4 sigma", part3()),
        ("2024 inputs: 2.2 sigma on quoted errors, 1.6 sigma with the pole-mass ambiguity, on either top mass", part4()),
        ("exact agreement would need m_t = 173.7 GeV", part5()),
        ("in c the same 1.6 sigma; sqrt(e) excluded by over a hundred sigma", part6()),
        ("tampers: matching dominates running; comp 123's budget set identified", part7()),
        ("PDG 2026 inputs move lambda_SM(M_*) by less than 0.0003 and the sigma by less than 0.3", part8()),
        ("run back down, e^-1/4 returns v = 247 GeV from m_h, 0.4 percent high; the old 0.01 percent is withdrawn", part9()),
    ]
    print()
    print("SUMMARY")
    for label, ok in results:
        print(f"  [{'PASS' if ok else 'FAIL'}] {label}")
    print()
    print("  CONCLUSION: the running is converged at two loops and the honest, current-input")
    print("  statement of the Higgs self-coupling test is lambda_SM(M_*) = 0.0930 +- 0.0008 on the")
    print("  2026 Review of Particle Physics (0.0932 +- 0.0008 on the 2024 edition) against")
    print("  e^-1/4 = 0.09197: consistent at 1.4 sigma, 1.2 percent, with the error dominated by the")
    print("  top mass and its pole-mass ambiguity.  The 1.0 sigma the book printed through v32.136")
    print("  rested on an anchor the literature's own formulas do not produce.")
    print("  ALL CHECKS PASSED" if all(ok for _, ok in results) else "  SOME CHECKS FAILED")


if __name__ == "__main__":
    main()
