"""
verify_P069.py — P69: F₄ Cubic Norm and the Jordan Exponential Identity

Central claim (Theorem 2.1):
    ln N_J(exp_J(X_sec)) = Tr(X_sec) = μ₀

where X_sec = diag(π, π², 4π³),  N_J(diag(a,b,c)) = abc,  exp_J(diag(a,b,c)) = diag(eᵃ,eᵇ,eᶜ).

Additional checks:
  - μ₀ = 4π³ + π² + π  (definition)
  - p₃(X_sec) = π·π²·4π³ = 4π⁶  (cubic norm of X_sec itself)
  - μ₁ = (16π³/5 + 3π²/4 + 2π/3)
  - MU = μ₁/μ₀
  - E_B/E_c = 4π²/(1+π)  (spectral ratio, paper claims ≈ 9.53)
"""

import mpmath

mpmath.mp.dps = 50

PASS = FAIL = 0
N = 0

def report(name, ok, claimed, actual):
    global PASS, FAIL, N
    N += 1
    ok = bool(ok)
    PASS += ok
    FAIL += (not ok)
    status = "PASS" if ok else "FAIL"
    print(f"  [{status}] {N:>2}. {name}")
    if not ok:
        print(f"          claimed={claimed}")
        print(f"          actual ={actual}")

pi = mpmath.pi

mu0 = 4*pi**3 + pi**2 + pi
mu1 = mpmath.mpf('16')*pi**3/5 + 3*pi**2/4 + 2*pi/3
MU  = mu1/mu0

# --- Check 1: μ₀ value ---
mu0_expected = mpmath.mpf('137.036')
ok1 = abs(mu0 - mu0_expected) < mpmath.mpf('0.001')
report("P69-1: μ₀ = 4π³+π²+π ≈ 137.036", ok1,
       "≈137.036", f"={float(mu0):.6f}")

# --- Check 2: The main cubic norm identity ---
# X_sec = diag(π, π², 4π³)
a, b, c = pi, pi**2, 4*pi**3

# Jordan exponential
exp_a, exp_b, exp_c = mpmath.exp(a), mpmath.exp(b), mpmath.exp(c)

# Cubic norm of exp_J(X_sec) = exp_a * exp_b * exp_c
N_expJ = exp_a * exp_b * exp_c

# log of cubic norm
ln_N_expJ = mpmath.log(N_expJ)

# Should equal Tr(X_sec) = μ₀
diff = abs(ln_N_expJ - mu0)
ok2 = (diff < mpmath.mpf('1e-40'))
report("P69-2: ln N_J(exp_J(X_sec)) = μ₀  (exact algebraic identity)",
       ok2,
       f"μ₀={float(mu0):.10f}",
       f"ln N_J(exp_J(X_sec))={float(ln_N_expJ):.10f}, diff={float(diff):.2e}")

# --- Check 3: p₃(X_sec) = 4π⁶ ---
p3_Xsec = a * b * c   # = π · π² · 4π³ = 4π⁶
p3_expected = 4*pi**6
diff3 = abs(p3_Xsec - p3_expected)
ok3 = (diff3 < mpmath.mpf('1e-40'))
report("P69-3: N_J(X_sec) = π·π²·4π³ = 4π⁶", ok3,
       f"4π⁶={float(p3_expected):.6f}", f"product={float(p3_Xsec):.6f}")

# --- Check 4: MU = μ₁/μ₀ ≈ 0.79334 ---
MU_expected = mpmath.mpf('0.793342')
ok4 = abs(MU - MU_expected) < mpmath.mpf('1e-5')
report("P69-4: MU = μ₁/μ₀ ≈ 0.79334", ok4,
       f"≈{float(MU_expected):.6f}", f"={float(MU):.6f}")

# --- Check 5: spectral ratio E_B/E_c = 4π²/(1+π) ≈ 9.53 ---
E_e = pi
E_b = pi**2
E_B = 4*pi**3
E_c = E_e + E_b
ratio = E_B / E_c   # = 4π³/(π+π²) = 4π²/(1+π)
ratio_formula = 4*pi**2 / (1+pi)
ratio_expected = mpmath.mpf('9.53')
ok5_a = abs(ratio - ratio_formula) < mpmath.mpf('1e-40')
ok5_b = abs(ratio - ratio_expected) < mpmath.mpf('0.01')
ok5 = ok5_a and ok5_b
report("P69-5: E_B/E_c = 4π²/(1+π) ≈ 9.53", ok5,
       f"≈{float(ratio_expected)}", f"={float(ratio):.4f}")

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