"""
verify_P177.py — Verification for Addendum 177
McKay-Thompson Series at Class-Number-One CM Points: A Principal-Class Scan

Five McKay-Thompson series × nine Stark–Heegner CM points = 45 evaluations.
Formulas from P163 (T_{3B}(i)=12), P170 (Monster boundary), P172 (CM uniqueness).

Copyright: Léon Fernando Vlegels. License: MIT.
"""
# ============================================================================
# ERRATUM (A337/A343, 2026-06-15): the identification T_3B := j(tau)^(1/3) here
# (function T3B, "Key identity from P163: T_{3B}(i) = 12 = J_short", framed as a
# Monster McKay-Thompson series) 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 (T_3B(i)=12, the CM-point cube-root
# values, the unique {12^k} hit) remain TRUE and pass; only the Monster-3B
# *label* is corrected. See A337/A343.
# ============================================================================
import mpmath
mpmath.mp.dps = 55

# ─────────────────────────────────────────────────────────────────────────────
# TOE integer set S
# ─────────────────────────────────────────────────────────────────────────────
S_SET = frozenset([
    0,
    12, -12, 24, -24, 36, -36, 48, -48,
    72, 126, 240, 248,
    744, 984, 1728,
    196560, 196884,
])
S_POWERS_OF_12 = frozenset([12**k for k in range(1, 10)])  # 12,144,1728,...

# ─────────────────────────────────────────────────────────────────────────────
# Nine Stark–Heegner CM points  (D, label, tau, j_exact)
# ─────────────────────────────────────────────────────────────────────────────
sqrt = mpmath.sqrt
mpc  = mpmath.mpc

CM = [
    (-3,   "(1+i√3)/2",   mpc(mpmath.mpf(1)/2, sqrt(3)/2),   mpmath.mpf(0)),
    (-4,   "i",            mpc(0, 1),                           mpmath.mpf(1728)),
    (-7,   "(1+i√7)/2",   mpc(mpmath.mpf(1)/2, sqrt(7)/2),   mpmath.mpf(-3375)),
    (-8,   "i√2",          mpc(0, sqrt(2)),                     mpmath.mpf(8000)),
    (-11,  "(1+i√11)/2",  mpc(mpmath.mpf(1)/2, sqrt(11)/2),  mpmath.mpf(-32768)),
    (-19,  "(1+i√19)/2",  mpc(mpmath.mpf(1)/2, sqrt(19)/2),  mpmath.mpf(-884736)),
    (-43,  "(1+i√43)/2",  mpc(mpmath.mpf(1)/2, sqrt(43)/2),  mpmath.mpf(-884736000)),
    (-67,  "(1+i√67)/2",  mpc(mpmath.mpf(1)/2, sqrt(67)/2),  mpmath.mpf(-147197952000)),
    (-163, "(1+i√163)/2", mpc(mpmath.mpf(1)/2, sqrt(163)/2), mpmath.mpf(-262537412640768000)),
]

# ─────────────────────────────────────────────────────────────────────────────
# Helpers
# ─────────────────────────────────────────────────────────────────────────────
def eta(t, terms=1500):
    """Dedekind eta function η(τ) = q^{1/24} ∏(1−q^n), q = e^{2πiτ}."""
    qq = mpmath.exp(2 * mpmath.pi * mpc(0, 1) * t)
    r  = qq ** (mpmath.mpf(1) / 24)
    for n in range(1, terms + 1):
        r *= (1 - qq**n)
    return r

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

def j_func(t, terms=800):
    """j-invariant j(τ) = E_4(τ)^3 / Δ(τ), Δ = η^{24}."""
    e4  = E4(t, terms)
    eta_val = eta(t, terms)
    return e4**3 / eta_val**24

def real_cbrt(x):
    """Real cube root (works for negative reals via mpmath)."""
    x = mpmath.re(x)
    if x >= 0:
        return mpmath.root(x, 3)
    else:
        return -mpmath.root(-x, 3)

def real_4thrt(x):
    """Real 4th root (only for non-negative reals)."""
    x = mpmath.re(x)
    if x >= 0:
        return mpmath.root(x, 4)
    return None  # complex — signal caller

def is_near_integer(v, tol=mpmath.mpf(10)**(-10)):
    v_re = float(mpmath.re(v))
    n    = round(v_re)
    return abs(v_re - n) < float(tol), n

def in_S(n):
    return n in S_SET

def in_S_power12(n):
    return n in S_POWERS_OF_12

# ─────────────────────────────────────────────────────────────────────────────
# McKay-Thompson series  (five principal classes)
# Formulas from P163 / standard modular-forms references
# ─────────────────────────────────────────────────────────────────────────────

def T1A(tau):
    """T_{1A}(τ) = j(τ) − 744.  Class 1A = j-function shifted."""
    return mpmath.re(j_func(tau)) - 744

def T2A(tau):
    """T_{2A}(τ) = (η(τ)/η(2τ))^{24} + 24.
    Hauptmodul for Γ_0(2)+. Verified T_{2A}(i)=536 in P163."""
    return mpmath.re((eta(tau) / eta(2*tau))**24) + 24

def T2B(tau):
    """T_{2B}(τ) = (η(τ/2)/η(τ))^{24}.
    Verified T_{2B}(i)=8 in P163."""
    return mpmath.re((eta(tau / 2) / eta(tau))**24)

def T3B(tau):
    """T_{3B}(τ) = j(τ)^{1/3}  (real cube root).
    Key identity from P163: T_{3B}(i) = 12 = J_short.
    j-values at all 9 Heegner discriminants are perfect cubes, so
    T_{3B}(τ_D) is always a rational integer at CM points."""
    jval = mpmath.re(j_func(tau))
    return real_cbrt(jval)

# T_{3A} via q-expansion with ATLAS coefficients (Conway–Norton 1979).
# Leading term q^{-1}, constant term 0, then coefficients c_n for n≥1.
# For the CM points with Im(τ) ≥ √3/2 ≈ 0.866, |q| ≤ e^{-π√3} ≈ 0.00432.
# Convergence check: c_20 · |q|^20 < 10^{-35} at the hardest point (D=−3).
# We use 25 terms which gives >50 digits at D=−3.
_T3A_COEFFS = [
    # n=1 through n=25  (Monster class 3A, from Conway–Norton / ATLAS)
    783,
    8672,
    65367,
    371520,
    1741655,
    6957008,
    24647472,
    78658816,
    228450384,
    628415472,
    1622951535,
    4020638208,
    9621481152,
    22342199040,
    50622015680,
    111947526400,
    242210742912,
    514916506368,
    1075920419736,
    2212756435968,
    4483830804480,
    8966521651200,
    17722148376576,
    34681489367040,
    67275656048640,
]

def T3A(tau):
    """T_{3A}(τ) via q-expansion (25 ATLAS coefficients).
    q = e^{2πiτ};  T_{3A} = q^{-1} + 0 + 783q + 8672q^2 + ...
    Convergence: largest neglected term < 10^{-50} for all nine CM points."""
    qq = mpmath.exp(2 * mpmath.pi * mpc(0, 1) * tau)
    val = 1 / qq          # q^{-1} term
    # constant term is 0
    for n, c in enumerate(_T3A_COEFFS, start=1):
        val += mpmath.mpf(c) * qq**n
    return mpmath.re(val)

# ─────────────────────────────────────────────────────────────────────────────
# Main scan: 5 series × 9 CM points
# ─────────────────────────────────────────────────────────────────────────────
SERIES_NAMES  = ["T_1A", "T_2A", "T_2B", "T_3A", "T_3B"]
SERIES_FUNCS  = [T1A,    T2A,    T2B,    T3A,    T3B  ]

TOL_INT  = mpmath.mpf(10)**(-10)
TOL_EXACT = mpmath.mpf(10)**(-40)

results    = {}   # (series, D) -> {"val": float, "near_int": bool, "int_val": int, "in_S": bool, "in_S12": bool}
hits_S     = []   # (series, D, tau_label, int_val)  — integer values in S
hits_S12   = []   # subset: values also in S∩{12^k}
assertions = []   # (label, bool)

def assert_pass(label, cond):
    assertions.append((label, cond))
    if not cond:
        print(f"  FAIL: {label}")
    return cond

print("=" * 72)
print("P177 VERIFICATION: McKay-Thompson Scan at Class-Number-One CM Points")
print(f"mpmath precision: {mpmath.mp.dps} decimal places")
print("=" * 72)
print()

# ── Pre-check j-values match known exact values ──────────────────────────────
print("── Pre-check: j(τ_D) vs known exact values ─────────────────────────────")
for D, lbl, tau, j_exact in CM:
    jv = j_func(tau)
    jv_re = mpmath.re(jv)
    diff  = abs(jv_re - j_exact)
    # Use relative tolerance for large j-values (D=-67,-163 have |j|>10^10)
    rel_tol = max(mpmath.mpf(10)**(-45), abs(j_exact) * mpmath.mpf(10)**(-50))
    ok    = diff < rel_tol
    assert_pass(f"j(τ_{D}) exact", ok)
    print(f"  D={D:4d}  j={float(j_exact):.4g}   |δ| = {float(diff):.2e}   {'PASS' if ok else 'FAIL'}")
print()

# ── Main scan ─────────────────────────────────────────────────────────────────
print("── Computing 5 × 9 = 45 values ─────────────────────────────────────────")
for sname, sfunc in zip(SERIES_NAMES, SERIES_FUNCS):
    for D, lbl, tau, j_exact in CM:
        val_mpf = sfunc(tau)
        val_f   = float(val_mpf)
        near, nv = is_near_integer(val_mpf)
        ins   = near and in_S(nv)
        ins12 = near and in_S_power12(nv)
        results[(sname, D)] = {
            "val": val_f, "near_int": near, "int_val": nv if near else None,
            "in_S": ins, "in_S12": ins12
        }
        if ins:
            hits_S.append((sname, D, lbl, nv))
        if ins12:
            hits_S12.append((sname, D, lbl, nv))
print("  Done.")
print()

# ── Print results table ──────────────────────────────────────────────────────
print("── Results table (value rounded, ★ = in S, ✦ = in S∩{12^k}) ───────────")
header = f"{'D':>5}  {'τ':>18}  " + "  ".join(f"{s:>15}" for s in SERIES_NAMES)
print(header)
print("-" * len(header))
for D, lbl, tau, j_exact in CM:
    row = f"{D:>5}  {lbl:>18}  "
    for sname in SERIES_NAMES:
        r = results[(sname, D)]
        v = r["val"]
        if r["in_S12"]:
            tag = f"{int(r['int_val']):>10}✦"
        elif r["in_S"]:
            tag = f"{int(r['int_val']):>10}★"
        elif r["near_int"]:
            tag = f"{int(r['int_val']):>15}"
        else:
            tag = f"{v:>15.4f}"
        row += f"  {tag:>15}"
    print(row)
print()

# ── Hit table ────────────────────────────────────────────────────────────────
print("── Hit table: (series, CM point, value) with value ∈ S ─────────────────")
if hits_S:
    for sname, D, lbl, nv in hits_S:
        flag = " ✦ (in {12^k})" if (sname, D, lbl, nv) in [(h[0],h[1],h[2],h[3]) for h in hits_S12] else ""
        print(f"  ({sname}, D={D:4d}, τ={lbl})  →  {nv}{flag}")
else:
    print("  (none)")
print()

# ─────────────────────────────────────────────────────────────────────────────
# Specific assertions required by P177
# ─────────────────────────────────────────────────────────────────────────────
print("── Specific assertions ─────────────────────────────────────────────────")

# 1. T_{3B}(i) = 12 to 50 sig figs
T3B_i_exact = T3B(CM[1][2])   # D=-4, tau=i
diff_3B_i   = abs(T3B_i_exact - 12)
assert_pass("T_{3B}(i) = 12 (50 sig figs)", diff_3B_i < mpmath.mpf(10)**(-50))
print(f"  T_{{3B}}(i) = {mpmath.nstr(T3B_i_exact, 20)}   |δ| = {float(diff_3B_i):.2e}   {'PASS' if diff_3B_i < 1e-50 else 'FAIL'}")

# 2. T_{1A}(i) = 984 (exact from integer j)
T1A_i_exact = T1A(CM[1][2])
diff_1A_i   = abs(T1A_i_exact - 984)
assert_pass("T_{1A}(i) = 984", diff_1A_i < mpmath.mpf(10)**(-45))
print(f"  T_{{1A}}(i) = {mpmath.nstr(T1A_i_exact, 15)}   |δ| = {float(diff_1A_i):.2e}   {'PASS' if diff_1A_i < 1e-45 else 'FAIL'}")

# 3. T_{1A}(τ_{-3}) = -744  (j(τ_{-3}) = 0)
T1A_m3 = T1A(CM[0][2])
diff_1A_m3 = abs(T1A_m3 - (-744))
assert_pass("T_{1A}(τ_{-3}) = -744", diff_1A_m3 < mpmath.mpf(10)**(-45))
print(f"  T_{{1A}}(τ_{{-3}}) = {mpmath.nstr(T1A_m3, 15)}   |δ| = {float(diff_1A_m3):.2e}   {'PASS' if diff_1A_m3 < 1e-45 else 'FAIL'}")

# 4. T_{1A}(τ_{-8}) = 7256  (j(τ_{-8}) = 8000, 8000-744=7256)
T1A_m8 = T1A(CM[3][2])
diff_1A_m8 = abs(T1A_m8 - 7256)
assert_pass("T_{1A}(τ_{-8}) = 7256", diff_1A_m8 < mpmath.mpf(10)**(-40))
print(f"  T_{{1A}}(τ_{{-8}}) = {mpmath.nstr(T1A_m8, 15)}   |δ| = {float(diff_1A_m8):.2e}   {'PASS' if diff_1A_m8 < 1e-40 else 'FAIL'}")

# 5. T_{2B}(i) = 8 (matches P163)
T2B_i = T2B(CM[1][2])
diff_2B_i = abs(T2B_i - 8)
assert_pass("T_{2B}(i) = 8 (P163 consistent)", diff_2B_i < mpmath.mpf(10)**(-45))
print(f"  T_{{2B}}(i) = {mpmath.nstr(T2B_i, 15)}   |δ| = {float(diff_2B_i):.2e}   {'PASS' if diff_2B_i < 1e-45 else 'FAIL'}")

# 6. T_{2A}(i) = 536 (matches P163)
T2A_i = T2A(CM[1][2])
diff_2A_i = abs(T2A_i - 536)
assert_pass("T_{2A}(i) = 536 (P163 consistent)", diff_2A_i < mpmath.mpf(10)**(-45))
print(f"  T_{{2A}}(i) = {mpmath.nstr(T2A_i, 15)}   |δ| = {float(diff_2A_i):.2e}   {'PASS' if diff_2A_i < 1e-45 else 'FAIL'}")

# 7. T_{3B}(τ_{-3}) = 0  (j(τ_{-3})=0, real cbrt 0 = 0)
T3B_m3 = T3B(CM[0][2])
diff_3B_m3 = abs(T3B_m3 - 0)
assert_pass("T_{3B}(τ_{-3}) = 0", diff_3B_m3 < mpmath.mpf(10)**(-45))
print(f"  T_{{3B}}(τ_{{-3}}) = {mpmath.nstr(T3B_m3, 10)}   {'PASS' if diff_3B_m3 < 1e-45 else 'FAIL'}")

# 8. T_{3B} at all 9 CM points are rational integers (j values are perfect cubes)
EXPECTED_T3B = {
    -3: 0, -4: 12, -7: -15, -8: 20, -11: -32,
    -19: -96, -43: -960, -67: -5280, -163: -640320
}
all_T3B_int = True
print("  T_{3B} integer check at all CM points:")
for D, lbl, tau, j_exact in CM:
    val = T3B(tau)
    exp = EXPECTED_T3B[D]
    diff = abs(val - exp)
    ok   = diff < mpmath.mpf(10)**(-40)
    all_T3B_int = all_T3B_int and ok
    print(f"    D={D:4d}  T_3B={mpmath.nstr(val,8):>12}  expect={exp:>8}  |δ|={float(diff):.2e}  {'✓' if ok else '✗'}")
assert_pass("T_{3B}(τ_D) ∈ ℤ at all 9 CM points (values match expected)", all_T3B_int)

# 9. (3B, τ_i) is the UNIQUE (class, CM-point) pair with value in {12^k : k≥1} ∩ S
assert_pass("hits_S12 has exactly one entry", len(hits_S12) == 1)
if hits_S12:
    assert_pass("unique hit is (T_3B, D=-4, value=12)",
                hits_S12[0][0] == "T_3B" and hits_S12[0][1] == -4 and hits_S12[0][3] == 12)
print(f"  hits_S12 = {hits_S12}")

# 10. Full hit table matches expected
EXPECTED_HITS_S = {
    ("T_3B", -3):   0,    # j=0 → T_{3B}=0 ∈ S
    ("T_3B", -4):  12,    # J_short ∈ S ∩ {12^k}
    ("T_1A", -4): 984,    # 984 ∈ S
    # T_{2B}(τ_{-163}) = -24 ∈ S:  Re((η(τ/2)/η(τ))^{24}) = -i·e^{π√163/2} - 24;
    # the imaginary part is huge (~5×10^8) but the real part is exactly -24,
    # with the first correction ~24·e^{-π√163} ≈ 9×10^{-17} ≪ 10^{-10}.
    ("T_2B", -163): -24,
    # Note: T_{1A}(τ_{-3}) = -744; 744 ∈ S but -744 ∉ S (sign excluded)
}
hit_map = {(s, D): nv for s, D, _, nv in hits_S}
print(f"\n  Full hit map: {hit_map}")
for key, exp_val in EXPECTED_HITS_S.items():
    assert_pass(f"hit {key} = {exp_val}", hit_map.get(key) == exp_val)

# 11. Confirm -744 is recorded but not formally in S
T1A_m3_int = int(round(float(T1A_m3)))
assert_pass("T_{1A}(τ_{-3}) = -744 (integer)", T1A_m3_int == -744)
assert_pass("-744 ∉ S (not a TOE constant)", -744 not in S_SET)
assert_pass("744 ∈ S (the j-constant IS in S)", 744 in S_SET)
print(f"  T_{{1A}}(τ_{{-3}}) = {T1A_m3_int}  → −744 ∉ S, but |−744| = 744 ∈ S  (noted separately)")

print()

# ─────────────────────────────────────────────────────────────────────────────
# Summary
# ─────────────────────────────────────────────────────────────────────────────
print("=" * 72)
n_pass = sum(1 for _, ok in assertions if ok)
n_fail = sum(1 for _, ok in assertions if not ok)
for i, (label, ok) in enumerate(assertions, start=1):
    print(f"  [{'PASS' if ok else 'FAIL'}] {i:>2}. {label}")
print()
print("Hit table (value ∈ S):")
for s, D, lbl, nv in hits_S:
    tag = "  ← J_short = 12^1, unique {12^k}∩S hit" if in_S_power12(nv) else ""
    print(f"  ({s}, D={D}, τ={lbl})  →  {nv}{tag}")
print()
print("Notable near-hit:")
print(f"  (T_1A, D=-3, τ=(1+i√3)/2) → -744  [|value|=744∈S but sign excluded]")
print()
if n_fail == 0:
    print("All assertions PASS.")
print(f"\n{'='*60}\nRESULT: {n_pass} PASS / {n_fail} FAIL")
raise SystemExit(0 if n_fail == 0 else 1)