"""
verify_P163.py — Verification script for Addendum 163
McKay-Thompson Series T_g(i) for Structurally Significant Monster Classes
L. F. Vlegels, May 2026. Copyright: Léon Fernando Vlegels. License: MIT.
"""
# ============================================================================
# ERRATUM (A337/A343, 2026-06-15): the identification T_3B := j(tau)^(1/3) here
# is RETRACTED. j^(1/3) = E4/eta^8 is the E8/G2 cube-root function; its value 12
# at tau=i is E8/G2 ROOT GEOMETRY (= J_short), NOT a Monster 3B McKay-Thompson
# value (canonical 3B = (eta/eta3)^12 + 12 = 535.59 at i). The arithmetic checks
# below (j(i)^(1/3) = 12, sections 5/8) remain TRUE and pass; only the
# Monster-3B *label* is corrected. See A337/A343.
# ============================================================================
import mpmath
mpmath.mp.dps = 50

tau = mpmath.mpc(0, 1)   # τ = i
q   = mpmath.exp(2 * mpmath.pi * mpmath.mpc(0, 1) * tau)   # q = e^{-2π}

# ──────────────────────────────────────────────────────────
# Helper: Dedekind eta at τ (product form)
# ──────────────────────────────────────────────────────────
def eta(t, terms=1000):
    qq = mpmath.exp(2 * mpmath.pi * mpmath.mpc(0, 1) * t)
    r  = qq ** (mpmath.mpf(1) / 24)
    for n in range(1, terms + 1):
        r *= (1 - qq**n)
    return r

# ──────────────────────────────────────────────────────────
# Helper: Eisenstein series E4 at τ
# ──────────────────────────────────────────────────────────
def E4(t, terms=500):
    qq = mpmath.exp(2 * mpmath.pi * mpmath.mpc(0, 1) * t)
    return mpmath.mpf(1) + 240 * sum(n**3 * qq**n / (1 - qq**n)
                                      for n in range(1, terms + 1))

# ──────────────────────────────────────────────────────────
# Precompute eta values at τ = i, 2i, i/2, 3i, 4i
# ──────────────────────────────────────────────────────────
eta_i    = eta(tau)
eta_2i   = eta(2 * tau)
eta_hi   = eta(tau / 2)
eta_3i   = eta(3 * tau)
eta_4i   = eta(4 * tau)

q_re = mpmath.re(q)   # real part; Im(q) ≈ 0 at τ = i

TOL_EXACT = mpmath.mpf(10) ** (-40)   # tolerance for claimed equalities
TOL_CLOSE = mpmath.mpf(10) ** (-6)    # tolerance for convergent series checks

PASS = FAIL = 0
N = 0

def record(label, ok):
    global PASS, FAIL, N
    ok = bool(ok)
    N += 1
    PASS += ok; FAIL += (not ok)
    print(f"  [{'PASS' if ok else 'FAIL'}] {N:>2}. {label}")
    return ok

def check(label, value, expected, tol=TOL_EXACT):
    diff = abs(value - expected)
    ok   = diff < tol
    record(label, ok)
    return ok, diff, "PASS" if ok else "FAIL"

print("=" * 70)
print("P163 VERIFICATION: McKay-Thompson T_g(i) for Monster Conjugacy Classes")
print("mpmath precision:", mpmath.mp.dps, "decimal places")
print("=" * 70)
print()

# ──────────────────────────────────────────────────────────
# § 0 — Basic CM values
# ──────────────────────────────────────────────────────────
print("── § 0  Basic CM values ──────────────────────────────────────────────")

j_i = mpmath.re(E4(tau)**3 / eta_i**24)
ok, diff, st = check("j(i) = 1728", j_i, 1728)
print(f"       j(i) = {mpmath.nstr(j_i, 15)}   |delta| = {float(diff):.3e}")

ok, diff, st = check("j(i)^{1/3} = 12", j_i**(mpmath.mpf(1)/3), 12)
print(f"       j(i)^{{1/3}} = {mpmath.nstr(j_i**(mpmath.mpf(1)/3), 15)}   |delta| = {float(diff):.3e}")
print()

# ──────────────────────────────────────────────────────────
# § 1 — T_1A(i) = J(i) = j(i) − 744 = 984
# ──────────────────────────────────────────────────────────
print("── § 1  T_{1A}(i) = J(i) = j(i) − 744 ──────────────────────────────")
T1A = j_i - 744
ok, diff, st = check("T_1A(i) = 984", T1A, 984)
print(f"       T_{{1A}}(i) = {mpmath.nstr(T1A, 15)}   (expect 984)")
print()

# ──────────────────────────────────────────────────────────
# § 2 — T_2A(i) = (η(i)/η(2i))^{24} + 24 = 536
# ──────────────────────────────────────────────────────────
print("── § 2  T_{2A}(i) — Hauptmodul for Γ₀(2)⁺ ──────────────────────────")
raw2A = mpmath.re((eta_i / eta_2i)**24)
ok, diff, st = check("(eta_i/eta_2i)^24 = 512 = 2^9", raw2A, 512)
print(f"       (η(i)/η(2i))^{{24}} = {mpmath.nstr(raw2A, 20)}   |delta| = {float(diff):.3e}")

T2A = raw2A + 24
ok, diff, st = check("T_2A(i) = 536", T2A, 536)
print(f"       T_{{2A}}(i) = (η/η₂)^24 + 24 = {mpmath.nstr(T2A, 20)}")

# CM identity: η(i)/η(2i) = 2^{3/8}
ratio_2A = mpmath.re(eta_i / eta_2i)
ok, diff, st = check("eta(i)/eta(2i) = 2^{3/8}", ratio_2A, mpmath.power(2, mpmath.mpf(3)/8))
print(f"       η(i)/η(2i) = 2^{{3/8}} = {mpmath.nstr(mpmath.power(2, mpmath.mpf(3)/8), 15)}")
print()

# ──────────────────────────────────────────────────────────
# § 3 — T_2B(i) = (η(i/2)/η(i))^{24} = 8
# ──────────────────────────────────────────────────────────
print("── § 3  T_{2B}(i) — via (η(i/2)/η(i))^{24} ─────────────────────────")
T2B = mpmath.re((eta_hi / eta_i)**24)
ok, diff, st = check("T_2B(i) = 8 = 2^3", T2B, 8)
print(f"       (η(i/2)/η(i))^{{24}} = {mpmath.nstr(T2B, 20)}   |delta| = {float(diff):.3e}")

# CM identity: η(i/2)/η(i) = 2^{1/8}
ratio_2B = mpmath.re(eta_hi / eta_i)
ok, diff, st = check("eta(i/2)/eta(i) = 2^{1/8}", ratio_2B, mpmath.power(2, mpmath.mpf(1)/8))
print(f"       η(i/2)/η(i) = 2^{{1/8}} = {mpmath.nstr(mpmath.power(2, mpmath.mpf(1)/8), 15)}")
print()

# ──────────────────────────────────────────────────────────
# § 4 — T_3A(i) — q-series evaluation (not an integer)
# ──────────────────────────────────────────────────────────
print("── § 4  T_{3A}(i) — q-expansion (7 terms) ───────────────────────────")
# Coefficients from Conway-Norton tables (class 3A)
coeffs_3A = {0: 0, 1: 783, 2: 8672, 3: 65367, 4: 371520, 5: 1741655, 6: 6957008}
T3A = 1 / q_re  # q^{-1} term
for k, c in coeffs_3A.items():
    T3A += c * q_re**k
print(f"  T_{{3A}}(i) ≈ {mpmath.nstr(T3A, 20)}")
print(f"  Not an integer (closest integer: {int(round(float(T3A.real if hasattr(T3A,'real') else T3A)))})")
is_int = abs(T3A - round(float(T3A))) > 0.01
record(f"T_3A non-integer: {'confirmed' if is_int else 'UNEXPECTED'}", is_int)
print()

# ──────────────────────────────────────────────────────────
# § 5 — T_3B(i) = j(i)^{1/3} = 12 = J_short
# ──────────────────────────────────────────────────────────
print("── § 5  T_{3B}(i) = j(i)^{1/3} = J_short ───────────────────────────")
T3B = j_i**(mpmath.mpf(1)/3)
ok, diff, st = check("T_3B(i) = 12 = J_short", T3B, 12)
print(f"       T_{{3B}}(i) = j(i)^{{1/3}} = {mpmath.nstr(T3B, 20)}   |delta| = {float(diff):.3e}")
print(f"       T_{{3B}}(i) = J_short = 12  (exact TOE constant)")
print()

# ──────────────────────────────────────────────────────────
# § 6 — T_4A(i) = j(i)^{1/4} = 12^{3/4} (algebraic irrational)
# ──────────────────────────────────────────────────────────
print("── § 6  T_{4A}(i) = j(i)^{1/4} = 12^{3/4} ──────────────────────────")
T4A = j_i**(mpmath.mpf(1)/4)
T4A_analytic = mpmath.power(12, mpmath.mpf(3)/4)
ok, diff, st = check("j(i)^{1/4} = 12^{3/4}", T4A, T4A_analytic)
print(f"       j(i)^{{1/4}} = {mpmath.nstr(T4A, 20)}")
print(f"  12^{{3/4}}   = {mpmath.nstr(T4A_analytic, 20)}")
print(f"  Not a TOE constant (irrational)")
print()

# ──────────────────────────────────────────────────────────
# § 7 — T_6A(i) = j(i)^{1/6} = √12 = 2√3 (algebraic irrational)
# ──────────────────────────────────────────────────────────
print("── § 7  T_{6A}(i) = j(i)^{1/6} = √12 ───────────────────────────────")
T6A = j_i**(mpmath.mpf(1)/6)
sqrt12 = mpmath.sqrt(12)
ok, diff, st = check("j(i)^{1/6} = sqrt(12)", T6A, sqrt12)
print(f"       j(i)^{{1/6}} = {mpmath.nstr(T6A, 20)}")
print(f"  √12 = 2√3   = {mpmath.nstr(sqrt12, 20)}")
print(f"  T_{{6A}}(i) = √J_short = √12  (algebraic, not the TOE constant itself)")
print()

# ──────────────────────────────────────────────────────────
# § 8 — Uniqueness: only 3B gives J_short exactly
# ──────────────────────────────────────────────────────────
print("── § 8  T_{3B} uniqueness: j(i)^{p/q} = J_short^{3p/q} ─────────────")
J_short = mpmath.mpf(12)
for p, q_exp, name in [(1,1,"1A"), (1,3,"3B"), (1,4,"4A"), (1,6,"6A")]:
    val = j_i**(mpmath.mpf(p)/q_exp)
    is_J = abs(val - J_short) < TOL_EXACT
    eq_str = f"= J_short" if is_J else f"≠ J_short"
    print(f"  j(i)^{{1/{q_exp}}} = {mpmath.nstr(val, 12):>22}  {eq_str}  (T_{{{name}}} pattern)")
    if name == "3B":
        ok_u, diff_u, st_u = check(f"T_3B uniqueness", val, J_short)
        record(f"j^{{1/3}} = J_short uniqueness", ok_u)

print(f"  Only j(i)^{{1/3}} = J_short (3p/q ∈ ℤ iff q|3 and p/q = 1/3)")
print()

# ──────────────────────────────────────────────────────────
# § 9 — TOE constant cross-check
# ──────────────────────────────────────────────────────────
print("── § 9  TOE constant cross-check ────────────────────────────────────")
ALPHA_INV = 4*mpmath.pi**3 + mpmath.pi**2 + mpmath.pi
BREATH    = mpmath.pi * ALPHA_INV
phi       = (1 + mpmath.sqrt(5)) / 2
TOE = {"J_short": 12, "ALPHA_INV": float(ALPHA_INV), "BREATH": float(BREATH),
       "B_G2": 48, "B_F4": 36, "j_i": 1728}
for name, val in [("T_1A=984", 984), ("T_2A=536", 536), ("T_2B=8", 8)]:
    hits = [k for k,v in TOE.items() if abs(v - val) < 0.5]
    match = hits if hits else ["none"]
    print(f"  {name}: TOE matches = {match}")
print()

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