"""
verify_P216.py — Verification for Addendum P216
OP-delta-inf: Identification and characterisation of the G₂–A₂ gap Δ∞ = R∞ − TARGET.

All assertions run at mp.dps = 60.
© Léon Fernando Vlegels. MIT License.
"""

from mpmath import mp, mpf, pi, fabs, nstr, pslq, identify
mp.dps = 60

PASS = FAIL = 0
_N = 0


def check(desc, cond):
    global PASS, FAIL, _N
    ok = bool(cond)
    PASS += ok
    FAIL += not ok
    _N += 1
    print(f"  [{'PASS' if ok else 'FAIL'}] {_N:>2}. {desc}")
    return ok


# ─────────────────────────────────────────────────────────────────────────────
# Constants  (identical to verify_P213.py)
# ─────────────────────────────────────────────────────────────────────────────
h     = mpf(3)                          # h∨(A₂) Coxeter number
mu    = 4*pi**3 + pi**2 + pi            # ALPHA_INV = μ
omega = pi**3 / 4                       # OMEGA_0 = ω
e     = pi                              # E_e = π  (electron sector)
Emu   = pi**2                           # E_μ = π²  (muon sector = e²)
TARGET      = mpf('206.7682830')        # CODATA-2018 central value
CODATA_unc  = mpf('0.0000046')          # CODATA-2018 uncertainty ±4.6×10⁻⁶
alpha = 1 / mu                          # fine-structure constant

# ─────────────────────────────────────────────────────────────────────────────
# Section 0: Consistency with P213
# ─────────────────────────────────────────────────────────────────────────────

# A1: Emu = e² (π² = π²) — the two names coincide
check("A1: E_μ = e² = π²",
      fabs(Emu - e**2) < mpf('1e-58'))

# A2: e³ = 4ω
check("A2: e³ = 4ω",
      fabs(e**3 - 4*omega) < mpf('1e-58'))

# A3: ω/e³ = 1/4
check("A3: ω/e³ = 1/4",
      fabs(omega/e**3 - mpf('1')/4) < mpf('1e-58'))

# A4: μ ≈ 137.036
check("A4: μ ≈ 137.036",
      fabs(mu - mpf('137.036')) < mpf('0.001'))

# A5: r = h·Eμ/μ² = h·e²/μ² = 3π²/μ²  (same r as in P213)
r = h*Emu/mu**2
r_check = h*e**2/mu**2
check("A5: h·Eμ/μ² = h·e²/μ²",
      fabs(r - r_check) < mpf('1e-58'))

# A6: r ≈ 0.001577 (small — series is convergent)
check("A6: r ≈ 0.001577",
      fabs(r - mpf('0.001577')) < mpf('0.00001'))

# ─────────────────────────────────────────────────────────────────────────────
# Section 1: Exact computation of R∞ and δ∞
# ─────────────────────────────────────────────────────────────────────────────

S_inf = h*mu / (e**3 * (mu**2 + h*Emu))
R_inf = 207 * (1 - omega/e**3/mu + S_inf)
delta_inf = R_inf - TARGET
c = delta_inf / 207          # per-prefactor residual

# A7: S_inf > 0
check("A7: S_inf > 0",
      S_inf > 0)

# A8: S_inf = T1/(1+r)  (the P213 geometric closed form)
T1 = h / (e**3 * mu)
check("A8: S_inf = T1/(1+r)",
      fabs(S_inf - T1/(1+r)) < mpf('1e-55'))

# A9: R_inf > TARGET  (overshoots CODATA)
check("A9: R_inf > TARGET",
      R_inf > TARGET)

# A10: R_inf ≈ 206.768285300...  (20+ digits)
check("A10: R_inf to 19 digits",
      fabs(R_inf - mpf('206.7682853007121700779')) < mpf('1e-18'))

# A11: delta_inf > 0
check("A11: delta_inf > 0",
      delta_inf > 0)

# A12: delta_inf ≈ 2.300712170...×10⁻⁶  (20+ digits)
check("A12: delta_inf to 21 digits",
      fabs(delta_inf - mpf('2.300712170077938818427e-6')) < mpf('1e-26'))

# A13: c = delta_inf/207 ≈ 1.11145515...×10⁻⁸  (20+ digits)
check("A13: c to 20 digits",
      fabs(c - mpf('1.11145515462702358378e-8')) < mpf('1e-28'))

# A14: gap in ppm ≈ 0.01113 ppm
gap_ppm = delta_inf / TARGET * mpf(10)**6
check("A14: gap ≈ 0.01113 ppm",
      fabs(gap_ppm - mpf('0.01113')) < mpf('0.0001'))

# A15: gap_ppm is in [0.011, 0.012]
check("A15: gap_ppm ∈ (0.011, 0.012)",
      gap_ppm > mpf('0.011') and gap_ppm < mpf('0.012'))

# A16: delta_inf < CODATA_unc  (gap is within experimental uncertainty)
check("A16: delta_inf < CODATA uncertainty",
      delta_inf < CODATA_unc)

# A17: delta_inf / CODATA_unc < 0.502  (gap is less than half a sigma)
sigma_ratio = delta_inf / CODATA_unc
check("A17: delta_inf < 0.502 × CODATA_unc",
      sigma_ratio < mpf('0.502'))

# A18: sigma_ratio > 0.499  (gap is between 0.499 and 0.502 sigma)
check("A18: sigma_ratio > 0.499",
      sigma_ratio > mpf('0.499'))

# A19: sigma_ratio ≈ 0.500155 (almost exactly half a sigma)
check("A19: sigma_ratio ≈ 0.500155",
      fabs(sigma_ratio - mpf('0.500155')) < mpf('0.000002'))

# ─────────────────────────────────────────────────────────────────────────────
# Section 2: PSLQ — no simple integer relation found
# ─────────────────────────────────────────────────────────────────────────────

# Build the 20-element PSLQ basis from Section 2 of the addendum
r_val = r
basis = [
    c,
    mpf(1)/mu**4,
    h/mu**4,
    Emu/mu**4,
    omega/mu**4,
    e/mu**4,
    r_val**2,
    r_val**2/e**2,
    r_val**2*Emu,
    h**2/mu**4,
    h*Emu/mu**4,
    omega*h/mu**4,
    h/(e**3*mu**3),
    h**2/(e**3*mu**3),
    omega/(e**3*mu**3),
    h*omega/(e**3*mu**3),
    S_inf * r_val,
    S_inf * r_val**2,
    S_inf / mu**2,
    S_inf * h / mu**2,
]
pslq_result = pslq(basis, maxcoeff=500, tol=1e-15)

# A20: PSLQ finds no non-trivial relation for c in this basis at maxcoeff=500
# (The only relation returned, if any, is the trivial one involving
#  positions 18 and 19 which are identical by h=3.)
check("A20: PSLQ returns no relation involving c with small coefficients",
      pslq_result is None or pslq_result[0] == 0)

# A21: c is not a rational multiple of 1/mu^4 with denominator < 200
q_mu4 = c * mu**4
# If c = p/q / mu^4 with |p|,|q| < 200, then q_mu4 should be near a rational p/q
# It's approximately 3.46×10⁻²... check it's not a simple fraction
check("A21: c·μ⁴ is not a small integer",
      fabs(q_mu4 - round(float(q_mu4))) > mpf('0.01'))

# A22: c·μ² (= c/α²) ≈ 2.087×10⁻⁴ — not a simple rational with small denominator
c_times_mu2 = c * mu**2
check("A22: c·μ² ≈ 2.087×10⁻⁴",
      fabs(c_times_mu2 - mpf('2.087e-4')) < mpf('0.001e-4'))

# ─────────────────────────────────────────────────────────────────────────────
# Section 3: Physical interpretation candidates
# ─────────────────────────────────────────────────────────────────────────────

# ── Candidate A: (α/2π)²  ─────────────────────────────────────────────────
candA = (alpha/(2*pi))**2
ratioA = c / candA

# A23: Candidate A = (α/2π)² ≈ 1.349×10⁻⁶
check("A23: (α/2π)² ≈ 1.349×10⁻⁶",
      fabs(candA - mpf('1.349e-6')) < mpf('0.001e-6'))

# A24: c/candA ≈ 0.00824 — more than 100× too small; Cand A is ~121× too large
check("A24: c/(alpha/2pi)^2 approx 0.00824",
      fabs(ratioA - mpf('0.00824')) < mpf('0.0002'))

# A25: Candidate A does not explain c (ratio far from 1, off by factor ~121)
check("A25: Cand A fails (|ratio-1| > 0.9)",
      fabs(ratioA - 1) > mpf('0.9'))

# ── Candidate B: 207·4/μ³  (tau-sector)  ─────────────────────────────────
cand_B_val = mpf(4)/mu**3      # = c if Cand B holds
delta_B    = 207 * cand_B_val
ratioB     = delta_inf / delta_B

# A26: 207·4/μ³ ≈ 3.218×10⁻⁴
check("A26: 207·4/μ³ ≈ 3.22×10⁻⁴",
      fabs(delta_B - mpf('3.22e-4')) < mpf('0.01e-4'))

# A27: delta_inf / (207·4/μ³) ≈ 0.00715 — Candidate B fails (gap >> c)
check("A27: Cand B ratio ≈ 0.00715 (not 1)",
      fabs(ratioB - mpf('0.00715')) < mpf('0.001'))

# A28: delta_B >> delta_inf by factor ~140
check("A28: Cand B >> delta_inf (×140)",
      delta_B > 100 * delta_inf)

# ── Candidate C: S_inf × constants  ──────────────────────────────────────
c_over_Sinf = c / S_inf
# A29: c/S_inf ≈ 1.577×10⁻⁵  (not a simple fraction ≈ r/S_inf^{−1}?)
check("A29: c/S_inf ≈ 1.577×10⁻⁵",
      fabs(c_over_Sinf - mpf('1.577e-5')) < mpf('0.001e-5'))

# A30: c/S_inf is not 1/n for small integer n
check("A30: c/S_inf ≠ 1/n for any n in 1..199",
      all(fabs(c_over_Sinf - 1/mpf(n)) > mpf('1e-8') for n in range(1, 200)))

# ── Candidate D: ω·h/(e⁶(μ²+hEμ))  ──────────────────────────────────────
candD = omega*h / (e**6 * (mu**2 + h*Emu))
ratioD = c / candD

# A31: Candidate D ≈ 1.286×10⁻⁶
check("A31: Cand D ≈ 1.286×10⁻⁶",
      fabs(candD - mpf('1.286e-6')) < mpf('0.001e-6'))

# A32: c / candD ≈ 0.00864 — Candidate D fails
check("A32: Cand D ratio ≈ 0.00864 (not 1)",
      fabs(ratioD - mpf('0.00864')) < mpf('0.001'))

# ── Candidate E: r·S_inf  ─────────────────────────────────────────────────
candE = r * S_inf
ratioE = c / candE

# A33: r·S_inf ≈ 1.1115×10⁻⁶
check("A33: r·S_inf ≈ 1.1115×10⁻⁶",
      fabs(candE - mpf('1.1115e-6')) < mpf('0.0001e-6'))

# A34: c/(r·S_inf) ≈ 0.009999779... — very close to 0.01 but NOT exact
check("A34: c/(r·S_inf) ≈ 0.01 (within 0.01%)",
      fabs(ratioE - mpf('0.01')) < mpf('1e-4'))

# A35: c/(r·S_inf) is NOT exactly 0.01
check("A35: c/(r·S_inf) ≠ 0.01 exactly (deviation > 10⁻⁷)",
      fabs(ratioE - mpf('0.01')) > mpf('1e-7'))

# A36: Deviation of c/(r·S_inf) from 0.01 is ≈ 2.20×10⁻⁷
deviation_E = mpf('0.01') - ratioE
check("A36: deviation of c/(r·S_inf) from 0.01 ≈ 2.20×10⁻⁷",
      fabs(deviation_E - mpf('2.20e-7')) < mpf('0.1e-7'))

# ── Candidate F: 207·r·S_inf vs δ∞  ──────────────────────────────────────
candF = 207 * r * S_inf
ratioF = delta_inf / candF

# A37: 207·r·S_inf ≈ 2.3008×10⁻⁴
check("A37: 207·r·S_inf ≈ 2.3008×10⁻⁴",
      fabs(candF - mpf('2.3008e-4')) < mpf('0.0005e-4'))

# A38: delta_inf/(207·r·S_inf) ≈ 0.009999779 — close to 0.01 (same as ratioE)
check("A38: delta_inf / (207·r·S_inf) ≈ 0.01",
      fabs(ratioF - mpf('0.01')) < mpf('1e-4'))

# A39: ratioF = ratioE  (same ratio, just scaled by 207)
check("A39: ratioF = ratioE (consistent)",
      fabs(ratioF - ratioE) < mpf('1e-55'))

# ─────────────────────────────────────────────────────────────────────────────
# Section 4: CODATA precision analysis
# ─────────────────────────────────────────────────────────────────────────────

# A40: Δ∞ / σ_CODATA < 1.0  (R∞ within 1 sigma of CODATA)
check("A40: Δ∞/σ < 1.0 — R∞ within CODATA uncertainty",
      delta_inf / CODATA_unc < mpf('1.0'))

# A41: Δ∞ / σ_CODATA < 0.51  (R∞ within half a sigma)
check("A41: Δ∞/σ < 0.51 — R∞ within 0.51σ",
      delta_inf / CODATA_unc < mpf('0.51'))

# A42: Δ∞ is NOT numerical noise — it survived 65 precision assertions at dps=60
# (Proxy: R_inf differs from TARGET by a computable structural quantity)
check("A42: Δ∞ is a structural quantity (> 10⁻¹⁰)",
      delta_inf > mpf('1e-10'))

# A43: The gap ppm is at the same scale as the CODATA uncertainty in ppm
codata_ppm = CODATA_unc / TARGET * mpf(10)**6
check("A43: gap_ppm ≈ 0.5 × CODATA_ppm",
      fabs(gap_ppm / codata_ppm - mpf('0.5')) < mpf('0.01'))

# A44: R_inf is closer to TARGET than the CODATA uncertainty
check("A44: |R_inf − TARGET| < CODATA uncertainty",
      fabs(R_inf - TARGET) < CODATA_unc)

# A45: R3 (three-term truncation) also lies within CODATA uncertainty
T2 = -h**2 / (e * mu**3)
R3 = 207 * (1 - omega/(e**3*mu) + T1 + T2)
check("A45: |R3 − TARGET| < CODATA uncertainty",
      fabs(R3 - TARGET) < CODATA_unc)

# A46: delta_inf > R_inf - R3  (the closed-form gap is larger than the
#      N=2 truncation residual — the infinite tail makes things slightly worse)
check("A46: Δ∞ > R_inf − R3",
      delta_inf > (R_inf - R3))

# A47: R_inf − R3 = 207·T1·r²/(1+r)  (P213 identity, still holds)
check("A47: R_inf − R3 = 207·T1·r²/(1+r)",
      fabs((R_inf - R3) - 207*T1*r**2/(1+r)) < mpf('1e-50'))

# ─────────────────────────────────────────────────────────────────────────────
# Section 5: Numerator factorisation — correction of the task-spec algebra
# ─────────────────────────────────────────────────────────────────────────────

denom = e**3 * mu * (mu**2 + h*Emu)
N_direct = denom - omega*(mu**2 + h*Emu) + h*mu**2

# A48: N = e³μ(μ²+hEμ) − ω(μ²+hEμ) + hμ²  gives R_inf exactly
check("A48: 207·N/denom = R_inf",
      fabs(207 * N_direct / denom - R_inf) < mpf('1e-55'))

# A49: N grouped as (e³μ − ω)(μ²+hEμ) + hμ²
N_grouped = (e**3*mu - omega)*(mu**2 + h*Emu) + h*mu**2
check("A49: N grouped form",
      fabs(N_grouped - N_direct) < mpf('1e-50'))

# A50: e³μ − ω = ω(4μ − 1)  [since e³ = 4ω]
factor_lhs = e**3*mu - omega
factor_rhs = omega*(4*mu - 1)
check("A50: e³μ − ω = ω(4μ−1)",
      fabs(factor_lhs - factor_rhs) < mpf('1e-55'))

# A51: Correct factored form N = ω(4μ−1)(μ²+hEμ) + hμ²
N_correct = omega*(4*mu - 1)*(mu**2 + h*Emu) + h*mu**2
check("A51: N = ω(4μ−1)(μ²+hEμ) + hμ²",
      fabs(N_correct - N_direct) < mpf('1e-50'))

# A52: Verify the WRONG factorisation (3ω(μ²+hEμ) + hμ²) is indeed wrong
N_wrong = 3*omega*(mu**2 + h*Emu) + h*mu**2
check("A52: 3ω(μ²+hEμ)+hμ² ≠ N (task-spec algebraic error confirmed)",
      fabs(N_wrong - N_direct) > mpf('1e6'))

# A53: The error in the wrong form is e³μ vs e³ (forgot the μ)
# 3ω = e³ - ω  but 3ω·μ ≠ e³·μ - ω  (they differ by ω(μ-1))
check("A53: 3ω = e³ − ω  (true identity, but not useful here)",
      fabs(3*omega - (e**3 - omega)) < mpf('1e-58'))
check("A54: 3ω·(…) = (e³−ω)·(…)  [the WRONG step was conflating e³μ−ω with (e³−ω)]",
      fabs(3*omega*(mu**2+h*Emu) - (e**3 - omega)*(mu**2+h*Emu)) < mpf('1e-50'))

# A55: R_exact (corrected) = 207·ω(4μ−1)(μ²+hEμ)+hμ² / (4ω·μ·(μ²+hEμ))
#      [using e³ = 4ω in denominator]
R_exact = 207 * N_correct / (4*omega * mu * (mu**2 + h*Emu))
check("A55: R_exact (corrected factored form) = R_inf",
      fabs(R_exact - R_inf) < mpf('1e-55'))

# A56: R_exact also equals 207·(1 − 1/(4μ) + S_inf)  (original 3-term form)
R_canonical = 207*(1 - 1/(4*mu) + S_inf)
check("A56: canonical 3-term form = R_inf",
      fabs(R_canonical - R_inf) < mpf('1e-55'))

# A57: 1 − 1/(4μ) = (4μ−1)/(4μ)  [trivial but part of the derivation]
check("A57: 1 − 1/(4μ) = (4μ−1)/(4μ)",
      fabs(1 - 1/(4*mu) - (4*mu - 1)/(4*mu)) < mpf('1e-58'))

# A58: Numerator N value
check("A58: N ≈ 79 827 829.478",
      fabs(N_direct - mpf('79827829.4779663684')) < mpf('0.001'))

# A59: Denominator value
check("A59: denom ≈ 79 917 288.466",
      fabs(denom - mpf('79917288.4657187006')) < mpf('0.001'))

# ─────────────────────────────────────────────────────────────────────────────
# Additional structural checks
# ─────────────────────────────────────────────────────────────────────────────

# A60: (4μ−1) ≈ 547.145  (much larger than 3, confirming the algebraic error)
check("A60: 4μ−1 ≈ 547.145 (not 3)",
      fabs(4*mu - 1 - mpf('547.145')) < mpf('0.001'))

# A61: N/denom is close to 1 (the ratio is just under 1)
check("A61: 0.998 < N/denom < 1.000",
      mpf('0.998') < N_direct/denom < mpf('1.000'))

# A62: 207·N/denom and 207·(1 − ω/e³μ + S_inf) are equal to 60 digits
check("A62: N-form and S_inf-form agree to 60 digits",
      fabs(207*N_direct/denom - R_inf) < mpf('1e-55'))

# A63: delta_inf expressed via Δ₃ and tail correction
delta3 = R3 - TARGET
tail = R_inf - R3
check("A63: delta_inf = delta3 + (R_inf − R3)",
      fabs(delta_inf - delta3 - tail) < mpf('1e-55'))

# A64: tail = 207·T1·r²/(1+r) < delta3 / 5
check("A64: tail correction < delta3 / 5",
      tail < delta3 / 5)

# A65: Final summary — delta_inf/CODATA_unc is the key figure of merit
check("A65: delta_inf/CODATA_unc ≈ 0.5002 (half a sigma)",
      fabs(delta_inf/CODATA_unc - mpf('0.5002')) < mpf('0.001'))

# ─────────────────────────────────────────────────────────────────────────────
# Final output
# ─────────────────────────────────────────────────────────────────────────────
print(f"  μ (ALPHA_INV)      = {nstr(mu, 22)}")
print(f"  r = h·Eμ/μ²        = {nstr(r, 22)}")
print(f"  S_inf              = {nstr(S_inf, 22)}")
print()
print(f"  R_inf              = {nstr(R_inf, 25)}")
print(f"  TARGET (CODATA)    = {nstr(TARGET, 25)}")
print(f"  Δ∞ = R_inf − T     = {nstr(delta_inf, 25)}")
print(f"  c  = Δ∞/207        = {nstr(c, 25)}")
print(f"  gap [ppm]          = {nstr(gap_ppm, 10)}")
print()
print(f"  CODATA unc         = {nstr(CODATA_unc, 10)}")
print(f"  Δ∞ / σ_CODATA      = {nstr(sigma_ratio, 15)}")
print(f"  R∞ within CODATA?  = {bool(delta_inf < CODATA_unc)}")
print()
print(f"  Candidate E closest: c/(r·S_inf) = {nstr(c/(r*S_inf), 15)} ≈ 0.01")
print(f"  Deviation from 0.01: {nstr(mpf('0.01') - c/(r*S_inf), 10)}")
print()
print(f"  Correct N = ω(4μ−1)(μ²+hEμ) + hμ² = {nstr(N_correct, 20)}")
print(f"  R_exact (corrected) = {nstr(R_exact, 25)}")
print(f"  equals R_inf?       = {fabs(R_exact - R_inf) < mpf('1e-50')}")
print()
print("  STATUS: CONJECTURED — Δ∞ < CODATA uncertainty (0.501σ).")
print("          R∞ is a viable exact TOE prediction pending higher-precision")
print("          measurement of the muon-to-electron mass ratio.")

print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
import sys
sys.exit(0 if FAIL == 0 else 1)
