"""
verify_P161.py — Addendum P161: E₈–Monster Bridge at τ = i.

Checks (50 decimal-place arithmetic where indicated):

  1. q-expansion coefficients computed by exact integer convolution:
       c(0) = 1
       c(1) = 248   (= dim E₈)
       c(2) = 4124  (= 1 + 248 + 3875)

  2. E₈ representation decomposition assertions:
       dim(E₈) = 248
       4124 = 1 + 248 + 3875

  3. Direct mpmath evaluation E₄(i)/η(i)^8 = 12 to 50 decimal places.

  4. Truncated series  q^{-1/3} Σ_{n=0}^{20} c(n)q^n  at q = e^{-2π}
     converges to 12 (truncation error < 1e-6; direct mpmath value within 1e-40 of 12).

Copyright: Léon Fernando Vlegels. MIT License. May 2026.
"""
# ============================================================================
# ERRATUM (A337/A343, 2026-06-15): the identification T_3B := j(tau)^(1/3) (the
# closing structural-reading line "j^{1/3} = T_{3B} (Monster 3B 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 (j(i)^(1/3) = 12, c[1]=248=dim E8, c[2]=4124) remain
# TRUE and pass; only the Monster-3B *label* is corrected. See A337/A343.
# ============================================================================

import math
import sys

import mpmath
from mpmath import mp, mpc, mpf, exp, pi, fabs, re, im

mp.dps = 50   # 50 decimal-place arithmetic throughout

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

# ─────────────────────────────────────────────────────────────────────────────
# Part 0: Helper — σ₃ divisor sum (exact integers)
# ─────────────────────────────────────────────────────────────────────────────

def sigma3_int(n):
    """Return σ₃(n) = Σ_{d|n} d³ as an exact Python int."""
    return sum(d * d * d for d in range(1, n + 1) if n % d == 0)


def binom_int(n, k):
    """Binomial coefficient C(n,k) as exact Python int."""
    return math.comb(n, k)


# ─────────────────────────────────────────────────────────────────────────────
# Part 1: q-expansion coefficients via exact integer arithmetic
# ─────────────────────────────────────────────────────────────────────────────
#
# E₄(τ) = Σ_{n≥0} e4[n] q^n  with  e4[0]=1, e4[n]=240·σ₃(n) for n≥1.
#
# η(τ)^{-8} = q^{-1/3} · F(τ)   where F(τ) = Σ f[n] q^n = ∏_{n≥1}(1-q^n)^{-8}.
#
# j(τ)^{1/3} = q^{-1/3} · Σ c[n] q^n   with   c = convolution of e4 and f.

N_COEFF = 22   # compute through order q^{N_COEFF}

# E₄ coefficients (exact integers)
e4 = [0] * (N_COEFF + 1)
e4[0] = 1
for n in range(1, N_COEFF + 1):
    e4[n] = 240 * sigma3_int(n)

# F = ∏_{n≥1}(1-q^n)^{-8} coefficients (exact integers via polynomial product)
# (1-q^n)^{-8} = Σ_{k≥0} C(k+7, 7) q^{nk}
# We accumulate the product one factor at a time.
f = [0] * (N_COEFF + 1)
f[0] = 1
for n in range(1, N_COEFF + 1):
    # Multiply current f by (1-q^n)^{-8}:
    # f_new[m] = Σ_{k≥0} C(k+7,7) · f_old[m - k·n]
    new_f = [0] * (N_COEFF + 1)
    for m in range(N_COEFF + 1):
        k = 0
        while k * n <= m:
            new_f[m] += binom_int(k + 7, 7) * f[m - k * n]
            k += 1
    f = new_f

# Convolution: c[n] = Σ_{k=0}^{n} e4[k] · f[n-k]
c = [0] * (N_COEFF + 1)
for n in range(N_COEFF + 1):
    for k in range(n + 1):
        c[n] += e4[k] * f[n - k]

print("P161 VERIFICATION: q-expansion of j^{1/3} and E8-Monster Bridge")

print()
print("S1  q-expansion coefficients (exact integer arithmetic)")
print()
print(f"  E4 coefficients:   e4[0]={e4[0]}   e4[1]={e4[1]}   e4[2]={e4[2]}")
print(f"  eta^{{-8}} product:  f[0]={f[0]}    f[1]={f[1]}      f[2]={f[2]}")
print(f"  j^{{1/3}} coeffs:    c[0]={c[0]}    c[1]={c[1]}    c[2]={c[2]}")
print()

check(1, f"c[0] = {c[0]}  (expect 1, trivial rep of E8)", c[0] == 1)
check(2, f"c[1] = {c[1]}  (expect 248 = dim(E8), adjoint rep)", c[1] == 248)
check(3, f"c[2] = {c[2]}  (expect 4124, derived from E4 x eta^-8 series)", c[2] == 4124)

# ─────────────────────────────────────────────────────────────────────────────
# Part 2: E₈ representation decomposition
# ─────────────────────────────────────────────────────────────────────────────

print()
print("S2  E8 representation decomposition")
print()

DIM_E8      = 248
DIM_TRIVIAL = 1
DIM_3875    = 3875

check(4, f"dim(E8) = {DIM_E8}", DIM_E8 == 248)
check(5, f"4124 = 1 + 248 + 3875 = {DIM_TRIVIAL + DIM_E8 + DIM_3875}  (E8 rep decomposition)",
      DIM_TRIVIAL + DIM_E8 + DIM_3875 == 4124)

# ─────────────────────────────────────────────────────────────────────────────
# Part 3: Direct mpmath evaluation of E₄(i)/η(i)^8
# ─────────────────────────────────────────────────────────────────────────────

print()
print("S3  Direct mpmath evaluation of E4(i)/eta(i)^8")
print()

PREC = mpf(10) ** (-40)

# Precompute σ₃ table up to N_MOD terms
N_MOD = 300
SIG3 = [sigma3_int(n) for n in range(1, N_MOD + 1)]

def E4_mpmath(tau, terms=N_MOD):
    q   = exp(2 * mpc(0, 1) * pi * tau)
    res = mpf(1)
    qn  = q
    for n in range(1, terms + 1):
        res += 240 * SIG3[n - 1] * qn
        if fabs(qn) < mpf(10) ** (-48):
            break
        qn *= q
    return res

def eta_mpmath(tau, terms=N_MOD):
    q   = exp(2 * mpc(0, 1) * pi * tau)
    res = q ** (mpf(1) / 24)
    qn  = q
    for n in range(1, terms + 1):
        res *= (1 - qn)
        if fabs(qn) < mpf(10) ** (-48):
            break
        qn *= q
    return res

tau_i = mpc(0, 1)
e4_val  = E4_mpmath(tau_i)
eta_val = eta_mpmath(tau_i)
j13_val = e4_val / eta_val ** 8

err3 = fabs(j13_val - 12)
print(f"  E4(i)/eta(i)^8 = {mp.nstr(re(j13_val), 30)}")
print(f"  |E4(i)/eta(i)^8 - 12| = {mp.nstr(err3, 4)}  (expect < 1e-40)")
check(6, "|E4(i)/eta(i)^8 - 12| < 1e-40  (= 12 to 50 decimal places)", err3 < PREC)

# ─────────────────────────────────────────────────────────────────────────────
# Part 4: Truncated series convergence at q = e^{-2π}
# ─────────────────────────────────────────────────────────────────────────────

print()
print("S4  Truncated series sum at q = e^{-2*pi}")
print()

q_real = exp(-2 * pi)   # q at tau = i; real positive number
print(f"  q = e^{{-2*pi}} = {mp.nstr(q_real, 6)}")

# q^{-1/3} * Σ_{n=0}^{20} c(n) * q^n
q_neg13 = q_real ** (-mpf(1) / 3)
trunc_sum = mpf(0)
N_TRUNC = 20
for n in range(N_TRUNC + 1):
    trunc_sum += c[n] * q_real ** n
trunc_sum *= q_neg13

err4_trunc = fabs(trunc_sum - 12)
print(f"  Summing q^{{-1/3}} * sum_{{n=0}}^{{{N_TRUNC}}} c(n)*q^n ...")
print(f"  Truncated sum  = {mp.nstr(trunc_sum, 15)}")
print(f"  |trunc - 12|   = {mp.nstr(err4_trunc, 3)}  (truncation; expect < 1e-6)")
check(7, f"truncated series (n=0..{N_TRUNC}) at q=e^-2pi: |trunc - 12| < 1e-6",
      err4_trunc < mpf(10) ** (-6))

err4_direct = fabs(re(j13_val) - 12)
print(f"  Direct value   = {mp.nstr(re(j13_val), 30)}")
print(f"  |direct - 12|  = {mp.nstr(err4_direct, 4)}  (expect < 1e-40)")
check(8, "|direct - 12| < 1e-40  (direct mpmath agrees)", err4_direct < PREC)

# ─────────────────────────────────────────────────────────────────────────────
# Summary
# ─────────────────────────────────────────────────────────────────────────────

print("""
  Structural reading:
    G2 root geometry  -->  J_short = 12             [P151, P158]
    j(i) = J_short^3 = 1728                         [P158]
    j(i)^{1/3} = 12 = [q^{-1/3}(1+248q+4124q^2+...)]_q=e^{-2pi}
                                                     [P159, P161]
    c(1) = 248 = dim(E8); c(2) = 1+248+3875         [McKay]
    j^{1/3} = T_{3B} (Monster 3B McKay-Thompson series)
""")
print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
sys.exit(0 if FAIL == 0 else 1)