"""
verify_P209.py — Verification for Addendum P209
Monad-electron-omega reformulation of candidate 8; PSLQ / second-term / 207 investigations.

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

from mpmath import mp, mpf, pi, nstr, log, exp, pslq, nint, sqrt, fabs
mp.dps = 60

# ─────────────────────────────────────────────────────────────────────────────
# Constants
# ─────────────────────────────────────────────────────────────────────────────
ALPHA_INV = 4*pi**3 + pi**2 + pi   # μ (monad)
OMEGA_0   = pi**3 / 4               # ω
E_e       = pi                       # e  (electron sector)
E_mu      = pi**2                    # muon sector
E_tau     = 4*pi**3                  # tau sector
TARGET    = mpf('206.7682830')       # CODATA muon/electron mass ratio

alpha     = 1 / ALPHA_INV           # fine-structure constant

PASS = FAIL = 0

def check(desc, cond):
    global PASS, FAIL
    ok = bool(cond); PASS += ok; FAIL += (not ok)
    print(f"  [{'PASS' if ok else 'FAIL'}] {PASS+FAIL:>2}. {desc}")
alpha_2pi = alpha / (2*pi)          # α/(2π) — first QED coupling

# Candidate 8 reformulated
g_alpha = OMEGA_0 / (E_e**3 * ALPHA_INV)
R_c8    = 207 * (1 - g_alpha)

# Angle B rank-1 formula
c2_rank1 = mpf(3) / (E_e**3 * ALPHA_INV)
R_rank1  = 207 * (1 - g_alpha + c2_rank1)

# ─────────────────────────────────────────────────────────────────────────────
# Section 1: Primitive identities
# ─────────────────────────────────────────────────────────────────────────────

# A1: ALPHA_INV is the correct fine-structure inverse at these digits
check("A1: μ sanity", abs(float(ALPHA_INV) - 137.036) < 0.001)

# A2: OMEGA_0 = π³/4 exactly
check("A2: ω = π³/4", abs(OMEGA_0 - E_e**3 / 4) < mpf('1e-58'))

# A3: ω/e³ = 1/4 exactly
check("A3: ω/e³ = 1/4", abs(OMEGA_0 / E_e**3 - mpf('1')/4) < mpf('1e-58'))

# A4: g_alpha = 1/(4μ)
check("A4: g_alpha = 1/(4μ)", abs(g_alpha - 1/(4*ALPHA_INV)) < mpf('1e-58'))

# A5: 4·g_alpha = 1/μ = α
check("A5: 4g = α", abs(4*g_alpha - alpha) < mpf('1e-58'))

# A6: E_tau = 4·e³
check("A6: E_tau = 4e³", abs(E_tau - 4*E_e**3) < mpf('1e-58'))

# A7: μ = E_tau + E_mu + E_e
check("A7: μ = E_tau+E_mu+E_e", abs(ALPHA_INV - (E_tau + E_mu + E_e)) < mpf('1e-58'))

# A8: ω = E_tau/16
check("A8: ω = E_tau/16", abs(OMEGA_0 - E_tau/16) < mpf('1e-58'))

# A9: α/(2π) = 1/(2π·μ)
check("A9: α/(2π) identity", abs(alpha_2pi - 1/(2*pi*ALPHA_INV)) < mpf('1e-58'))

# A10: c₂ = 6/(E_mu · 2π · μ) — alternative expression
c2_alt = 6 / (E_mu * 2 * pi * ALPHA_INV)
check("A10: c₂ = 6/(E_mu·2π·μ)", abs(c2_alt - c2_rank1) < mpf('1e-58'))

# ─────────────────────────────────────────────────────────────────────────────
# Section 2: Candidate 8 properties
# ─────────────────────────────────────────────────────────────────────────────

# B1: R_c8 is approximately 206.622
check("B1: R_c8 ≈ 206.622", abs(float(R_c8) - 206.622) < 0.01)

# B2: R_c8 gap from TARGET ~0.0706%
gap_c8 = abs(R_c8 - TARGET) / TARGET * 100
check("B2: R_c8 gap ≈ 0.0706%", abs(float(gap_c8) - 0.0706) < 0.001)

# B3: R_c8 < TARGET (undershoots)
check("B3: R_c8 < TARGET", R_c8 < TARGET)

# B4: R_c8 is closer to TARGET than 207 (improvement over tree-level)
gap_207 = abs(mpf(207) - TARGET) / TARGET * 100
check("B4: R_c8 improves on 207", gap_c8 < gap_207)

# B5: g_alpha < 1 (small coupling)
check("B5: g_alpha is small", g_alpha < mpf('0.01'))

# B6: g_alpha > 0
check("B6: g_alpha > 0", g_alpha > 0)

# B7: Δ = TARGET/207 - (1-g_alpha) > 0 (correction is positive)
Delta = TARGET / 207 - (1 - g_alpha)
check("B7: Δ > 0", Delta > 0)

# B8: Δ ≈ 7.05e-4
check("B8: Δ ≈ 7.05e-4", abs(float(Delta) - 7.05e-4) < 1e-5)

# B9: formula is purely in {ω, e, μ}: reformulation check
R_c8_check = 207 * (1 - OMEGA_0 / (E_e**3 * ALPHA_INV))
check("B9: reformulation consistency", abs(R_c8_check - R_c8) < mpf('1e-58'))

# B10: The correction 207·g_alpha ≈ 0.377
corr = 207 * g_alpha
check("B10: 207·g_alpha ≈ 0.377", abs(float(corr) - 0.377) < 0.01)

# ─────────────────────────────────────────────────────────────────────────────
# Section 3: Angle A — PSLQ null results (verify trivial identities found)
# ─────────────────────────────────────────────────────────────────────────────

# C1: Trivial basis identity: ω/e³ - e³/μ - ω/(e·μ) - ω/(e²·μ) = 0
trivial_1 = OMEGA_0/E_e**3 - E_e**3/ALPHA_INV - OMEGA_0/(E_e*ALPHA_INV) - OMEGA_0/(E_e**2*ALPHA_INV)
check("C1: trivial basis identity 1 vanishes", abs(trivial_1) < mpf('1e-55'))

# C2: Second trivial identity: 4ω = e³ (equivalently ω/e³ = 1/4)
trivial_2 = 4*OMEGA_0 - E_e**3
check("C2: 4ω = e³", abs(trivial_2) < mpf('1e-57'))

# C3: Third trivial identity: 1/μ² = 3ω/(e³μ²) + 4ω²/(e⁶μ²)
trivial_3 = 1/ALPHA_INV**2 - 3*OMEGA_0/(E_e**3*ALPHA_INV**2) - 4*OMEGA_0**2/(E_e**6*ALPHA_INV**2)
check("C3: third trivial identity vanishes", abs(trivial_3) < mpf('1e-60'))

# C4: Δ is NOT zero (i.e., PSLQ null result is meaningful — true gap exists)
check("C4: Δ is non-trivially large", abs(Delta) > mpf('1e-4'))

# C5: Δ is not expressible as ω/e³ (= 1/4) — sanity check
check("C5: Δ ≠ ω/e³", abs(Delta - OMEGA_0/E_e**3) > mpf('1e-2'))

# ─────────────────────────────────────────────────────────────────────────────
# Section 4: Angle B — second-term results
# ─────────────────────────────────────────────────────────────────────────────

# D1: c₂ = 3/(e³μ) > 0 (positive, pushes R toward TARGET)
check("D1: c₂ > 0", c2_rank1 > 0)

# D2: c₂ ≈ 7.06e-4
check("D2: c₂ ≈ 7.06e-4", abs(float(c2_rank1) - 7.06e-4) < 1e-5)

# D3: c₂ > g_alpha ... no, c₂ is smaller. Check c₂ < g_alpha:
# g_alpha = 1/(4μ) ≈ 1.824e-3, c₂ ≈ 7.06e-4. So c₂ < g_alpha.
check("D3: c₂ < g_alpha", c2_rank1 < g_alpha)

# D4: R_rank1 > R_c8 (rank-1 is closer to TARGET from below)
check("D4: R₁ > R_c8", R_rank1 > R_c8)

# D5: R_rank1 > TARGET (overshoots slightly)
check("D5: R₁ overshoots TARGET", R_rank1 > TARGET)

# D6: Gap of rank-1 formula < 0.0002% (well within 0.1%)
gap_rank1 = abs(R_rank1 - TARGET) / TARGET * 100
check("D6: gap < 0.0002%", float(gap_rank1) < 2e-4)

# D7: Gap of rank-1 formula is around 0.000112%
check("D7: gap ≈ 0.000112%", abs(float(gap_rank1) - 1.12e-4) < 1e-5)

# D8: Factor of improvement over R_c8
improvement = abs(R_c8 - TARGET) / abs(R_rank1 - TARGET)
check("D8: rank-1 improves R_c8 by >600x", float(improvement) > 600)

# D9: R_rank1 in combined form: 207·(1 + (3-ω)/(e³μ))
R_combined = 207 * (1 + (3 - OMEGA_0)/(E_e**3 * ALPHA_INV))
check("D9: combined form matches", abs(R_combined - R_rank1) < mpf('1e-50'))

# D10: c₂ = 6/(E_mu) · α/(2π)  [QED-structure connection]
check("D10: c₂ = 6/E_mu·α/(2π)", abs(6/E_mu * alpha_2pi - c2_rank1) < mpf('1e-55'))

# D11: The second-best formula (ω·e^{-4}·μ^{-1}) has gap > 0.012%
c2_rank2 = OMEGA_0 / (E_e**4 * ALPHA_INV)
R_rank2  = 207 * (1 - g_alpha + c2_rank2)
gap_rank2 = abs(R_rank2 - TARGET) / TARGET * 100
check("D11: rank-2 gap > 0.01%", float(gap_rank2) > 0.01)

# D12: Rank-1 formula is over 100x better than rank-2
check("D12: rank-1 >> rank-2", float(gap_rank2) / float(gap_rank1) > 100)

# D13: In the extended grid (n·ω^a·e^b·μ^c, n=1..6, a∈{0,1,2}, b∈{-6..1}, c∈{-2..1}),
# exactly one candidate sits within 0.001% of TARGET
candidates_within_001 = 0
best_in_grid = []
for n_int in range(1, 7):
    for a in [0, 1, 2]:
        for b in range(-6, 2):
            for c in [-2, -1, 0, 1]:
                c2_try = n_int * (OMEGA_0**a) * (E_e**b) * (ALPHA_INV**c)
                R_try = 207 * (1 - g_alpha + c2_try)
                gap_try = abs(R_try - TARGET) / TARGET * 100
                if gap_try < mpf('0.001'):
                    candidates_within_001 += 1
                    best_in_grid.append((gap_try, n_int, a, b, c))
check(f"D13: exactly 1 candidate within 0.001% (found {candidates_within_001}: {best_in_grid})", candidates_within_001 == 1)

# ─────────────────────────────────────────────────────────────────────────────
# Section 5: Angle C — 207 in {ω, e, μ}
# ─────────────────────────────────────────────────────────────────────────────

# E1: 207/μ ≈ 1.5105 (not 3/2 = 1.5 exactly)
ratio_207_mu = mpf(207) / ALPHA_INV
check("E1: 207/μ ≈ 1.5105", abs(float(ratio_207_mu) - 1.5105) < 0.001)

# E2: 3μ/2 ≠ 207 (best simple rational times μ misses by >0.5%)
err_3mu2 = abs(3*ALPHA_INV/2 - 207) / 207 * 100
check("E2: 3μ/2 ≠ 207", float(err_3mu2) > 0.5)

# E3: 19μ/(4π) is the best small-integer approximation, gap ~0.094%
approx_19 = 19 * ALPHA_INV / (4 * pi)
err_19 = abs(approx_19 - 207) / 207 * 100
check("E3: 19μ/(4π) within 0.1% of 207", float(err_19) < 0.1)

# E4: But 19μ/(4π) is too far for an ESTABLISHED formula
check("E4: 19μ/(4π) gap > 0.05% (not clean enough)", float(err_19) > 0.05)

# E5: μ + 70 ≈ 207 to 0.018%
err_mu70 = abs(ALPHA_INV + 70 - 207) / 207 * 100
check("E5: μ+70 within 0.02% of 207", float(err_mu70) < 0.02)

# E6: But again not a clean TOE expression
check("E6: μ+70 gap > 0.01% (70 is not a TOE primitive)", float(err_mu70) > 0.01)

# E7: 207 - round(μ) = 70 (confirming the residual)
check("E7: round(μ) = 137", int(nint(ALPHA_INV)) == 137)
check("E7b: 207 - 137 = 70", 207 - int(nint(ALPHA_INV)) == 70)

# E8: No combination e^a·μ^b·ω^c with |a|,|b|,|c|≤2 reproduces 207 within 0.1%
found_clean = False
for a_exp in range(-2, 3):
    for b_exp in range(-2, 3):
        for c_exp in range(-2, 3):
            if a_exp == 0 and b_exp == 0 and c_exp == 0:
                continue
            val = E_e**a_exp * ALPHA_INV**b_exp * OMEGA_0**c_exp
            if abs(float(val) - 207) / 207 < 0.001:
                found_clean = True
check("E8: no e^a·μ^b·ω^c combination reproduces 207 within 0.1%", not found_clean)

# E9: rank-1 formula expressed without explicit 207 requires it via the integer prefactor
# i.e., TARGET ≈ f(ω,e,μ)/g(ω,e,μ) only if 207 has a TOE form (which it doesn't)
# Verify that removing 207 and using μ directly doesn't work
R_mu_attempt = ALPHA_INV * (1 + (3 - OMEGA_0)/(E_e**3 * ALPHA_INV))
err_mu_attempt = abs(R_mu_attempt - TARGET) / TARGET * 100
check("E9: replacing 207→μ breaks the formula catastrophically", float(err_mu_attempt) > 30)

# ─────────────────────────────────────────────────────────────────────────────
# Section 6: Additional cross-checks
# ─────────────────────────────────────────────────────────────────────────────

# F1: α < 1/137 (fine structure is small)
check("F1: α < 1/136.9", float(alpha) < 1/136.9)
check("F1b: α > 1/137.1", float(alpha) > 1/137.1)

# F2: g_alpha = α/4
check("F2: g_alpha = α/4", abs(g_alpha - alpha/4) < mpf('1e-58'))

# F3: c₂ = 6α/(E_mu·2π) = 6/(E_mu·2π·μ) = 6/(π²·2π·μ) = 3/(π³μ) = 3/(e³μ) — consistency
c2_via_alpha = 6 * alpha / (E_mu * 2 * pi)
check("F3: c₂ = 6α/(E_mu·2π)", abs(c2_via_alpha - c2_rank1) < mpf('1e-55'))

# F4: The ratio c₂/g_alpha = 12/π³ (= 12/e³)
ratio_c2_g = c2_rank1 / g_alpha
check("F4: c₂/g_alpha = 12/π³", abs(ratio_c2_g - 12/E_e**3) < mpf('1e-55'))

# F5: 12/e³ = 12/π³ ≈ 0.3870
check("F5: 12/π³ ≈ 0.3870", abs(float(12/E_e**3) - 0.3870) < 0.001)

# F6: R_rank1 satisfies the 1 ppm threshold check (just barely above)
gap_ppm = abs(R_rank1 - TARGET) / TARGET * 1e6
check("F6: rank-1 gap > 1 ppm (OPEN status)", float(gap_ppm) > 1.0)
check("F6b: rank-1 gap < 1.2 ppm (remarkable near-miss)", float(gap_ppm) < 1.2)

# F7: R_c8 is within 1% of TARGET
check("F7: R_c8 within 1% of TARGET", float(abs(R_c8 - TARGET) / TARGET * 100) < 1.0)

# F8: Rank-1 is within 1 ppm of 1 ppm (within factor 2 of sub-ppm)
check("F8: rank-1 gap within factor 2 of 1 ppm", float(gap_ppm) < 2.0)

# F9: ω can be written as E_tau/16
check("F9: ω = E_tau/16", abs(OMEGA_0 - E_tau/16) < mpf('1e-58'))

# F10: The combined 3-term formula has correct sign structure
# 207·(1 - g_alpha + c₂) should overshoot TARGET slightly
check("F10: rank-1 formula overshoots TARGET", R_rank1 > TARGET)

# F11: TARGET is between R_c8 and R_rank1
check("F11: R_c8 < TARGET < R_rank1", R_c8 < TARGET < R_rank1)

# F12: c₂ ≈ Δ (the correction c₂ is very close to the needed gap Δ)
# Δ = TARGET/207 - (1 - g_alpha), c₂ = 3/(e³μ)
# They should be within 0.1% of each other
check("F12: c₂ ≈ Δ within 0.2%", abs(float(c2_rank1 - Delta) / float(Delta)) < 0.002)

# F13: c₂ - Δ ≈ 1.12e-6 (the residual after applying rank-1 correction)
residual = c2_rank1 - Delta
check("F13: |c₂ - Δ| < 2e-6", abs(float(residual)) < 2e-6)
check("F13b: c₂ > Δ (rank-1 overshoots)", float(residual) > 0)

# F14: Three-primitive structure: g_alpha uses exactly {ω, e, μ}
# g_alpha = ω / (e³ · μ) — no other primitives
g_recomputed = OMEGA_0 / (E_e**3 * ALPHA_INV)
check("F14: g_alpha recomputation", abs(g_recomputed - g_alpha) < mpf('1e-58'))

# F15: The QED connection: c₂ = 6/E_mu · α/(2π) suggests a connection to
# the muon mass sector (E_mu = π²) — verify the identity holds exactly
c2_from_qed = (6 / E_mu) * (1 / (2 * pi * ALPHA_INV))
check("F15: QED connection identity", abs(c2_from_qed - c2_rank1) < mpf('1e-57'))

print(f"\nKey numerical results (mp.dps={mp.dps}):")
print(f"  ALPHA_INV (μ) = {nstr(ALPHA_INV, 20)}")
print(f"  OMEGA_0   (ω) = {nstr(OMEGA_0, 20)}")
print(f"  E_e       (e) = {nstr(E_e, 20)}")
print(f"  g_alpha       = {nstr(g_alpha, 20)}")
print(f"  R_c8          = {nstr(R_c8, 20)}")
print(f"  R_rank1       = {nstr(R_rank1, 20)}")
print(f"  TARGET        = {nstr(TARGET, 20)}")
print(f"  gap_c8        = {nstr(gap_c8, 10)}%")
print(f"  gap_rank1     = {nstr(gap_rank1, 10)}% ({nstr(gap_rank1*1e4, 6)} ppm)")
print(f"  c₂_rank1      = {nstr(c2_rank1, 20)}")
print(f"  Δ             = {nstr(Delta, 20)}")
print(f"  improvement   = {nstr(improvement, 6)}x")
print(f"  19μ/(4π)      = {nstr(19*ALPHA_INV/(4*pi), 12)}, err={nstr(err_19, 6)}%")

print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
if FAIL:
    raise SystemExit(1)
