"""
verify_P206.py  —  Verification script for Addendum P206
"Mass Ratio from 4-Simplex Amplitude Ratios"

Systematically scans geometric quantities of the regular 4-simplex
(vertex inner products, geodesic arc lengths, face volumes, natural
coupling hypotheses) as candidates for the exponential coupling in
the sector-energy mass formula.

Author: Léon Fernando Vlegels
© Léon Fernando Vlegels. MIT License.
"""

from mpmath import mp, mpf, pi, exp, sqrt, acos, log, nstr, fabs

mp.dps = 60

# ── Constants ──────────────────────────────────────────────────────────────
ALPHA_INV  = 4*pi**3 + pi**2 + pi          # α⁻¹ ≈ 137.036
MU_0       = ALPHA_INV                      # μ₀ = ∫₀¹ ρ(x) dx
MU_1       = mpf(16)*pi**3/5 + 3*pi**2/4 + 2*pi/3  # μ₁ = ∫₀¹ x·ρ(x) dx
MU_A36     = MU_1 / MU_0                   # mass lever ≈ 0.7933
OMEGA_0    = pi**3 / 4                     # Ω₀ = π³/4

# Sector energies (Addendum 36, Proposition 2.1)
E_e        = pi                             # S¹/edge mode
E_mu       = pi**2                          # S³/boundary mode
E_tau      = 4*pi**3                        # B⁴/bulk mode
dE         = E_mu - E_e                    # π² − π ≈ 6.728

# CODATA 2018
CODATA     = mpf("206.7682830")

# Reference predictions from prior addenda
pred_sector = exp(MU_A36 * dE)             # A36 formula: 208.016
MU_fit      = log(mpf(207)) / dE           # coupling that gives exactly 207
pred_fit    = exp(MU_fit * dE)             # = 207.0 by construction

# ── 4-Simplex geometry (Addenda P202/P203) ─────────────────────────────────
theta4       = acos(mpf(-1) / 4)           # geodesic angle = arccos(−1/4) ≈ 1.8235 rad
chord        = sqrt(mpf(5) / 2)            # edge length √(5/2) in ℝ⁴
inner_prod   = mpf(-1) / 4                 # ⟨vᵢ, vⱼ⟩ for i≠j
gram_ev      = mpf(5) / 4                  # non-zero Gram eigenvalue (h*(A4)/h*(G2))

# Face/cell volumes with edge ℓ = √(5/2)
ell          = chord
V_face       = sqrt(3) / 4 * ell**2        # equilateral triangle area (type A₂)
V_cell       = sqrt(2) / 12 * ell**3       # regular tetrahedron volume (type A₃)

# Natural simplex coupling: MU_simplex = θ₄ / (E_μ − E_e)
MU_simplex   = theta4 / dE                 # ≈ 0.2710

# ── Helper ─────────────────────────────────────────────────────────────────
PASS = FAIL = 0
_N = 0

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

def rel_gap(x):
    return fabs(x - CODATA) / CODATA

def abs_gap(x):
    return fabs(x - CODATA)

# ── Group A: direct exponentials of geometric quantities ───────────────────
A1 = exp(theta4)                           # exp(arccos(−1/4))
A2 = exp(pi * theta4)                      # exp(π·arccos(−1/4))
A3 = A2                                    # same expression, labelled separately in paper
A4 = exp(theta4**2)                        # exp(θ₄²)
A5 = exp(theta4 * ALPHA_INV)              # exp(θ₄/α) — expected huge
A6 = exp(theta4 * ALPHA_INV / pi)         # exp(θ₄·α⁻¹/π)

# ── Group B: ratios with sector energies ───────────────────────────────────
B1 = exp(MU_simplex * dE)                  # = exp(θ₄) by definition
B2 = exp(gram_ev * dE)                     # exp((5/4)·(π²−π))
B3 = exp(gram_ev * pi * (pi - 1))          # exp((5/4)·π·(π−1))

# ── Group C: face / cell volumes ───────────────────────────────────────────
C4 = exp(V_cell)
C5 = exp(V_face)
C_ratio = V_cell / V_face

# ── Group D: natural coupling hypothesis ───────────────────────────────────
# Find continuous N such that exp(N·θ₄) = CODATA
N_exact = log(CODATA) / theta4             # ≈ 2.924 (non-integer)

# Evaluate at integer/rational/Lie-algebraic candidates
D_candidates = [
    ("N = h*(A1) = 2",          mpf(2)),
    ("N = h*(A2) = 3",          mpf(3)),
    ("N = h*(A3) = 4",          mpf(4)),
    ("N = h*(A4) = 5",          mpf(5)),
    ("N = 5/2",                 mpf(5) / 2),
    ("N = 8/3",                 mpf(8) / 3),
    ("N = 11/4",                mpf(11) / 4),
    ("N = π",                   pi),
    ("N = sqrt(10)",            sqrt(10)),
]

# ── Group E: combinations with Ω₀ and α⁻¹ ────────────────────────────────
E1 = exp(theta4 * OMEGA_0)                 # exp(θ₄·Ω₀) — expected huge
E2 = mpf(207) * (1 - theta4 / (4 * pi**2))  # 207·(1 − θ₄/(4π²))

# Additive series α⁻¹ + θ₄·Ω₀·k
E3 = {k: ALPHA_INV + theta4 * OMEGA_0 * k for k in range(-3, 4) if k != 0}

# ── OUTPUT ─────────────────────────────────────────────────────────────────
SEP = "=" * 72

print(SEP)
print("verify_P206.py  —  Mass Ratio from 4-Simplex Amplitude Ratios")
print(f"mpmath dps = {mp.dps}")
print(SEP)

print("\nFUNDAMENTAL CONSTANTS")
print(f"  α⁻¹         = {nstr(ALPHA_INV, 15)}")
print(f"  MU_A36      = μ₁/μ₀ = {nstr(MU_A36, 15)}")
print(f"  Ω₀          = π³/4  = {nstr(OMEGA_0, 15)}")
print(f"  CODATA      = m_μ/m_e = {nstr(CODATA, 12)}")

print("\n4-SIMPLEX GEOMETRY")
print(f"  θ₄ = arccos(−1/4) = {nstr(theta4, 15)} rad")
print(f"  θ₄ in degrees     = {nstr(theta4 * 180 / pi, 12)}°")
print(f"  chord  = √(5/2)   = {nstr(chord, 15)}")
print(f"  V_face = (√3/4)·(5/2)           = {nstr(V_face, 15)}")
print(f"  V_cell = (√2/12)·(5/2)^(3/2)    = {nstr(V_cell, 15)}")
print(f"  V_cell / V_face                  = {nstr(C_ratio, 15)}")
print(f"  MU_simplex = θ₄/(π²−π)          = {nstr(MU_simplex, 15)}")
print(f"  MU_A36 / MU_simplex              = {nstr(MU_A36 / MU_simplex, 15)}")
print(f"  N_exact = ln(CODATA)/θ₄          = {nstr(N_exact, 15)}")

print("\nREFERENCE PREDICTIONS")
print(f"  Sector formula (A36): {nstr(pred_sector, 12)}  (gap {nstr(rel_gap(pred_sector)*100,5)}%)")
print(f"  MU_fit = ln(207)/ΔE:  {nstr(pred_fit, 12)}  (gap {nstr(rel_gap(pred_fit)*100,5)}%)")

header = f"\n{'Candidate':<40}  {'Value':>18}  {'|Δ|':>12}  {'rel%':>9}"
rule   = "-" * 84
print(header)
print(rule)

all_candidates = []

def report(label, val):
    g  = abs_gap(val)
    rg = rel_gap(val)
    all_candidates.append((rg, label, val))
    print(f"  {label:<38}  {nstr(val, 10):>18}  {nstr(g, 6):>12}  {nstr(rg*100, 5):>9}%")

print("\nGROUP A — direct exponentials of geometric quantities")
report("A1  exp(θ₄)",                A1)
report("A2  exp(π·θ₄)",              A2)
report("A4  exp(θ₄²)",               A4)
# A5 and A6 are astronomically large; report informally
print(f"  {'A5  exp(θ₄·α⁻¹)':<38}  {'>> 10^100':>18}  (skip)")
print(f"  {'A6  exp(θ₄·α⁻¹/π)':<38}  {nstr(A6, 6):>18}  (informational: >> CODATA)")

print("\nGROUP B — ratios with sector energies")
report("B1  exp(MU_simplex·ΔE) = exp(θ₄)",  B1)
report("B2  exp((5/4)·(π²−π))",             B2)
report("B3  exp((5/4)·π·(π−1))",            B3)

print("\nGROUP C — face/cell volumes")
print(f"  {'C1  V_face = (√3/4)·(5/2)':<38}  {nstr(V_face,10):>18}  (informational)")
print(f"  {'C2  V_cell = (√2/12)·(5/2)^3/2':<38}  {nstr(V_cell,10):>18}  (informational)")
print(f"  {'C3  V_cell/V_face':<38}  {nstr(C_ratio,10):>18}  (informational)")
report("C4  exp(V_cell)",             C4)
report("C5  exp(V_face)",             C5)

print("\nGROUP D — natural coupling N·θ₄ hypothesis")
print(f"  Continuous N_exact = ln(CODATA)/θ₄ ≈ {nstr(N_exact, 10)}  (non-integer)")
for label, N_val in D_candidates:
    v = exp(N_val * theta4)
    report(f"D  exp(({label})·θ₄)", v)

print("\nGROUP E — combinations with Ω₀ and α⁻¹")
print(f"  {'E1  exp(θ₄·Ω₀)':<38}  {'>> 10^6':>18}  (informational: {nstr(E1,6)})")
report("E2  207·(1 − θ₄/(4π²))",    E2)
for k, v in sorted(E3.items()):
    print(f"  {'E3k={:+d}  α⁻¹+θ₄·Ω₀·{:+d}'.format(k,k):<38}  {nstr(v,10):>18}  (additive, informational)")

print()
print(rule)
print("TOP 5 CLOSEST TO CODATA (by relative gap):")
all_candidates.sort()
for rank, (rg, label, val) in enumerate(all_candidates[:5], 1):
    flag = " ← BEST" if rank == 1 else ""
    print(f"  #{rank}  {label:<44}  {nstr(val,10):>12}  rel={nstr(rg*100,5)}%{flag}")

print()
print(rule)
print("RATIO ANALYSIS  (is MU_A36 a simplex-integer multiple of MU_simplex?)")
print(f"  MU_A36 / MU_simplex  = {nstr(MU_A36/MU_simplex, 15)}")
print(f"  ln(CODATA) / θ₄      = {nstr(N_exact, 15)}")
print(f"  ln(207) / θ₄         = {nstr(log(mpf(207))/theta4, 15)}")
print(f"  All three ≈ 2.924, closest Lie integer = h*(A2) = 3")
print(f"  exp(3·θ₄)            = {nstr(exp(3*theta4), 12)}  (14.9% above CODATA)")

print()
print(rule)
print("COULD θ₄ EXPLAIN THE 0.60% MUON RESIDUAL?")
ratio_residual = pred_sector / CODATA
k_needed = theta4 / (1 - CODATA / pred_sector)
print(f"  208.016 / CODATA                 = {nstr(ratio_residual, 12)}")
print(f"  ln(208.016/CODATA)               = {nstr(log(ratio_residual), 12)}")
print(f"  k s.t. 208.016·(1 − θ₄/k) = CODATA: k = {nstr(k_needed, 12)}")
print(f"  k / 4π²                          = {nstr(k_needed/(4*pi**2), 12)}")
print(f"  Not a clean simplex number.")

print()
print(rule)
print("VERDICT")
best_rg = all_candidates[0][0]
best_label = all_candidates[0][1]
best_val   = all_candidates[0][2]
threshold_tight = mpf("0.0001")   # 0.01% → ESTABLISHED
threshold_loose = mpf("0.001")    # 0.1%  → PROMISING
if best_rg < threshold_tight:
    verdict = "ESTABLISHED (< 0.01% gap)"
elif best_rg < threshold_loose:
    verdict = "PROMISING (< 0.1% gap)"
else:
    verdict = "OPEN"
print(f"  Best simplex-only result: {best_label}")
print(f"  Value: {nstr(best_val, 10)}  relative gap: {nstr(best_rg*100, 5)}%")
print(f"  Threshold for ESTABLISHED: < 0.01%;  for OPEN: ≥ 0.1%")
print(f"  STATUS: {verdict}")
check("verdict: best simplex-only candidate gap >= 0.1% — STATUS OPEN "
      "(simplex geometry alone does not recover m_mu/m_e)",
      best_rg >= threshold_loose)
print()
print("  The 4-simplex geometry alone does not recover m_μ/m_e to within")
print("  the required precision.  The mass lever MU = μ₁/μ₀ encodes")
print("  independent physical content (depth-weighted density moment) that")
print("  is not determined by the simplex vertex structure.")

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