"""verify_P169.py — Verification for Addendum 169: The Leech Lattice Lambda_24.

Checks:
  1.  196560 = 2^4 * 3^3 * 5 * 7 * 13  (correct prime factorisation)
  2.  196560 / 240 = 819 = 3*273 = 9*91; 273 = 3*91; 91 = 7*13
  3.  Golay code dimension = 12 = J_short
  4.  Leech dimension = 24 = B_{G2}/2 = 48/2
  5.  24 = exponent in Delta(tau) = eta(tau)^24  (from P164)
  6.  3 * rank(E8) = 3 * 8 = 24  (Leech = three E8 shadows)
  7.  Theta series: [q^1] = 0 (no roots), [q^2] = 196560 (minimal vectors)
  8.  |Co_1| = 4157776806543360000  (exact integer, Conway group order)
  9.  |Co_1| = 2^21 * 3^9 * 5^4 * 7^2 * 11 * 13 * 23
 10.  196560 * J_short = 196560 * 12 = 2358720

Plus numerical verification:
     Theta_{Lambda_24}(i) = 1008 * eta(i)^24 = (j(i) - 720) * Delta(i)
     (two independent methods, relative error < 1e-50)

Copyright: Léon Fernando Vlegels. License: MIT. May 2026.
"""

from mpmath import mp, mpf, pi, exp, gamma, fabs, power, nstr, log10
import math
import sys

mp.dps = 55  # 55 digits → 50 sig-fig safety margin

# ── TOE constants ─────────────────────────────────────────────────────────────
ALPHA_INV   = 4*pi**3 + pi**2 + pi   # ≈ 137.036
J_short     = mpf(12)
j_i         = J_short**3             # = 1728
B_G2        = mpf(48)                # |Phi_{F4}| = B_{G2}  (P164)
B_F4        = mpf(36)                # proved P165

# Leech lattice invariants
LEECH_DIM   = 24            # dimension of Lambda_24
MIN_VECTORS = 196560        # number of minimal vectors (norm^2 = 4)
GOLAY_LEN   = 24            # length of the binary Golay code
GOLAY_DIM   = 12            # dimension of the binary Golay code
GOLAY_WMIN  = 8             # minimum Hamming weight

# E8 data (from P168)
RANK_E8     = 8
ROOTS_E8    = 240

# Co_1 (Conway group) order
CO1_ORDER   = 4157776806543360000

PASS = FAIL = 0
_N = 0

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

# ── 1. Prime factorisation of 196560 ─────────────────────────────────────────
print("=== 1. Prime factorisation: 196560 = 2^4 * 3^3 * 5 * 7 * 13 ===")

# Manual factorisation
def factorise(n):
    """Return dict {prime: exponent} for positive integer n."""
    factors = {}
    d = 2
    while d * d <= n:
        while n % d == 0:
            factors[d] = factors.get(d, 0) + 1
            n //= d
        d += 1
    if n > 1:
        factors[n] = factors.get(n, 0) + 1
    return factors

f = factorise(MIN_VECTORS)
print(f"  factorise(196560) = {f}")
check("196560 has prime factor 2 with exponent 4", f.get(2, 0) == 4)
check("196560 has prime factor 3 with exponent 3", f.get(3, 0) == 3)
check("196560 has prime factor 5 with exponent 1", f.get(5, 0) == 1)
check("196560 has prime factor 7 with exponent 1", f.get(7, 0) == 1)
check("196560 has prime factor 13 with exponent 1", f.get(13, 0) == 1)
check("196560 has exactly 5 distinct prime factors", len(f) == 5)

# Reconstruct and verify
reconstructed = 2**4 * 3**3 * 5 * 7 * 13
print(f"  2^4 * 3^3 * 5 * 7 * 13 = {reconstructed}")
check("2^4 * 3^3 * 5 * 7 * 13 = 196560", reconstructed == MIN_VECTORS)

# Primes 7 and 13 appear in Monster order
monster_primes = {2,3,5,7,11,13,17,19,23,29,31,41,47,59,71}
check("7 is a Monster prime", 7 in monster_primes)
check("13 is a Monster prime", 13 in monster_primes)

# ── 2. Ratio 196560 / 240 = 819 ───────────────────────────────────────────────
print("\n=== 2. Ratio 196560 / 240 = 819 = 3*273 = 9*91 = 9*7*13 ===")
ratio = MIN_VECTORS // ROOTS_E8
print(f"  196560 / 240 = {ratio}")
check("196560 / 240 = 819", ratio == 819)
check("819 = 3 * 273", ratio == 3 * 273)
check("273 = 3 * 91", 273 == 3 * 91)
check("91 = 7 * 13", 91 == 7 * 13)
check("819 = 9 * 91", ratio == 9 * 91)
check("819 = 9 * 7 * 13", ratio == 9 * 7 * 13)
check("196560 = 240 * 819", MIN_VECTORS == ROOTS_E8 * ratio)

f819 = factorise(819)
print(f"  factorise(819) = {f819}")
check("819 = 3^2 * 7 * 13", f819 == {3: 2, 7: 1, 13: 1})

# ── 3. Golay code dimension = 12 = J_short ───────────────────────────────────
print("\n=== 3. Golay code dimension = 12 = J_short ===")
print(f"  Golay code: [{GOLAY_LEN}, {GOLAY_DIM}, {GOLAY_WMIN}]")
check("Golay code dimension = 12", GOLAY_DIM == 12)
check("Golay code dimension = J_short", GOLAY_DIM == int(J_short))
check("Golay code length = 24 = Leech dim", GOLAY_LEN == LEECH_DIM)
check("Golay code minimum weight = 8", GOLAY_WMIN == 8)
check("Golay code: self-dual (dim = len/2)", GOLAY_DIM == GOLAY_LEN // 2)

# ── 4. Leech dimension = 24 = B_G2/2 ─────────────────────────────────────────
print("\n=== 4. Leech dimension = 24 = B_G2/2 = 48/2 ===")
bg2_half = int(B_G2) // 2
print(f"  B_G2 = {int(B_G2)};  B_G2/2 = {bg2_half}")
check("Leech dim = 24", LEECH_DIM == 24)
check("B_G2/2 = 24", bg2_half == 24)
check("Leech dim = B_G2/2", LEECH_DIM == bg2_half)
check("B_G2 = 48", int(B_G2) == 48)

# ── 5. 24 = exponent in Delta(tau) = eta(tau)^24 ─────────────────────────────
print("\n=== 5. 24 = exponent in Delta = eta^24 (from P164) ===")
ETA_EXPONENT = 24   # the exponent in Delta(tau) = eta(tau)^24
print(f"  Delta(tau) = eta(tau)^{ETA_EXPONENT}")
check("Eta exponent = 24", ETA_EXPONENT == 24)
check("Eta exponent = Leech dim", ETA_EXPONENT == LEECH_DIM)
check("Eta exponent = B_G2/2", ETA_EXPONENT == int(B_G2) // 2)

# Also: tau(2) = -24 (Ramanujan tau at 2)
tau_2 = -24   # Ramanujan tau function value
check("tau(2) = -24 = -Leech_dim", tau_2 == -LEECH_DIM)

# ── 6. 3 * rank(E8) = 24 ─────────────────────────────────────────────────────
print("\n=== 6. 3 * rank(E8) = 3 * 8 = 24 ===")
triple_e8 = 3 * RANK_E8
print(f"  3 * rank(E8) = 3 * {RANK_E8} = {triple_e8}")
check("3 * rank(E8) = 24", triple_e8 == 24)
check("3 * rank(E8) = Leech dim", triple_e8 == LEECH_DIM)
check("24 = 3 * 8 (explicit)", 24 == 3 * 8)

# The three-fold E8 structure: Leech in R^24 = R^8 + R^8 + R^8
print(f"  Leech in R^24 = R^8 ⊕ R^8 ⊕ R^8  (three E8 ambient spaces)")
check("3 * E8 root count = 3 * 240 = 720 (Niemeier E8^3 root count)", 3 * ROOTS_E8 == 720)

# ── 7. Theta series: [q^1]=0, [q^2]=196560 ───────────────────────────────────
print("\n=== 7. Theta series: [q^1]=0 (no roots), [q^2]=196560 ===")

def sigma3(n):
    return sum(d**3 for d in range(1, n+1) if n % d == 0)

# Fourier coefficients of E4^3 (up to q^4)
# E4 = 1 + 240*sigma3(1)*q + 240*sigma3(2)*q^2 + ...
e4_coeffs = [1] + [240 * sigma3(n) for n in range(1, 6)]

e8_coeffs = [0] * 6   # E4^2 = E8
for i in range(6):
    for j in range(6):
        if i + j < 6:
            e8_coeffs[i + j] += e4_coeffs[i] * e4_coeffs[j]

e43_coeffs = [0] * 6   # E4^3 = E4 * E8
for i in range(6):
    for j in range(6):
        if i + j < 6:
            e43_coeffs[i + j] += e4_coeffs[i] * e8_coeffs[j]

# Ramanujan tau function (first 5 values)
tau_vals = [0, 1, -24, 252, -1472, 4830]   # tau(0)=0, tau(1)=1, tau(2)=-24, ...

# Theta_{Lambda_24} = E4^3 - 720*Delta
leech_theta_coeffs = [e43_coeffs[n] - 720 * tau_vals[n] for n in range(5)]

print(f"  E4^3 coefficients [q^0..q^4]: {e43_coeffs[:5]}")
print(f"  Delta coefficients [q^0..q^4]: {tau_vals[:5]}")
print(f"  Theta_Leech [q^0..q^4]: {leech_theta_coeffs}")

check("[q^0] Theta_Leech = 1 (one zero vector)", leech_theta_coeffs[0] == 1)
check("[q^1] Theta_Leech = 0 (no roots = rootless)", leech_theta_coeffs[1] == 0)
check("[q^2] Theta_Leech = 196560 (minimal vectors)", leech_theta_coeffs[2] == 196560)
check("[q^3] Theta_Leech = 16773120 (second shell)", leech_theta_coeffs[3] == 16773120)
check("[q^4] Theta_Leech = 398034000 (third shell)", leech_theta_coeffs[4] == 398034000)

# The formula [q^2]: 179280 - 720*(-24) = 179280 + 17280 = 196560
q2_e43    = e43_coeffs[2]
q2_delta  = tau_vals[2]   # = -24
q2_leech  = q2_e43 - 720 * q2_delta
print(f"  [q^2] check: {q2_e43} - 720*({q2_delta}) = {q2_e43} + {-720*q2_delta} = {q2_leech}")
check("[q^2] from formula: 179280 - 720*(-24) = 196560", q2_leech == 196560)
check("[q^2] Theta_Leech = 196560 = MIN_VECTORS", q2_leech == MIN_VECTORS)

# ── 8. |Co_1| = 4157776806543360000 ──────────────────────────────────────────
print("\n=== 8. |Co_1| = 4157776806543360000 (Conway group order) ===")
co1_computed = 2**21 * 3**9 * 5**4 * 7**2 * 11 * 13 * 23
print(f"  |Co_1| [from factorisation] = {co1_computed}")
print(f"  |Co_1| [stated constant]    = {CO1_ORDER}")
check("|Co_1| matches stated value", co1_computed == CO1_ORDER)
check("|Co_1| = 4157776806543360000", CO1_ORDER == 4157776806543360000)
check("|Co_1| > 0", CO1_ORDER > 0)
check("|Co_0| = 2 * |Co_1|", 2 * CO1_ORDER == 2**22 * 3**9 * 5**4 * 7**2 * 11 * 13 * 23)

# ── 9. |Co_1| = 2^21 * 3^9 * 5^4 * 7^2 * 11 * 13 * 23 ───────────────────────
print("\n=== 9. Factorisation |Co_1| = 2^21 * 3^9 * 5^4 * 7^2 * 11 * 13 * 23 ===")
f_co1 = factorise(CO1_ORDER)
print(f"  factorise(|Co_1|) = {f_co1}")
check("|Co_1| factor 2^21", f_co1.get(2, 0) == 21)
check("|Co_1| factor 3^9",  f_co1.get(3, 0) == 9)
check("|Co_1| factor 5^4",  f_co1.get(5, 0) == 4)
check("|Co_1| factor 7^2",  f_co1.get(7, 0) == 2)
check("|Co_1| factor 11^1", f_co1.get(11, 0) == 1)
check("|Co_1| factor 13^1", f_co1.get(13, 0) == 1)
check("|Co_1| factor 23^1", f_co1.get(23, 0) == 1)
check("|Co_1| has exactly 7 distinct prime factors", len(f_co1) == 7)

# Verify all prime factors of Co_1 divide |Monster|
co1_primes = set(f_co1.keys())
monster_prime_factors = {2,3,5,7,11,13,23}   # those appearing in Co_1
check("All Co_1 prime factors are Monster primes",
      co1_primes.issubset(monster_primes))

# ── 10. 196560 * J_short = 2358720 ────────────────────────────────────────────
print("\n=== 10. 196560 * J_short = 196560 * 12 = 2358720 ===")
product = MIN_VECTORS * int(J_short)
print(f"  196560 * 12 = {product}")
check("196560 * 12 = 2358720", product == 2358720)
check("196560 * J_short = 2358720", product == 2358720)
check("2358720 / 12 = 196560", 2358720 // 12 == MIN_VECTORS)
check("2358720 = 196560 * J_short", 2358720 == MIN_VECTORS * int(J_short))

# Cross-check: 2358720 = 9801 * 240 + ...? Just confirm integer consistency.
check("2358720 divisible by 240", 2358720 % 240 == 0)
check("2358720 / 240 = 9828 = 819 * 12", 2358720 // 240 == 819 * 12)

# ── Numerical: Theta_{Lambda_24}(i) = 1008 * eta(i)^24 ───────────────────────
print("\n=== Numerical: Theta_Leech(i) = 1008 * eta(i)^24 [50 sig figs] ===")
q = exp(-2 * pi)   # q = e^{-2pi} at tau = i

# Method A: closed-form via eta(i) = Gamma(1/4) / (2 * pi^(3/4))
eta_i_cf = gamma(mpf('1') / 4) / (2 * power(pi, mpf('3') / 4))
Delta_i_cf = eta_i_cf**24
E4_i_cf = J_short * eta_i_cf**8        # E4(i) = 12 * eta(i)^8  (P168)
Theta_i_formula = (j_i - 720) * Delta_i_cf   # = 1008 * eta(i)^24

print(f"  eta(i) [closed form] = {nstr(eta_i_cf, 20)}")
print(f"  eta(i)^24 = Delta(i) = {nstr(Delta_i_cf, 20)}")
print(f"  E4(i) = 12*eta(i)^8  = {nstr(E4_i_cf, 20)}")
print(f"  j(i) - 720 = 1728 - 720 = {int(j_i - 720)}")
print(f"  Theta_Leech(i) [formula] = {nstr(Theta_i_formula, 20)}")

# Method B: via q-expansion of E4^3 - 720*Delta
N_TERMS = 30

# E4(i) via q-expansion
E4_i_qexp = mpf(1)
for n in range(1, N_TERMS + 1):
    E4_i_qexp += 240 * sigma3(n) * q**n

# Delta(i) via q-product: q * prod_{n>=1} (1-q^n)^24
Delta_i_qprod = q
for n in range(1, N_TERMS + 10):
    Delta_i_qprod *= (1 - q**n)**24
    if q**n < mpf('1e-65'):
        break

Theta_i_qexp = E4_i_qexp**3 - 720 * Delta_i_qprod
print(f"  Theta_Leech(i) [q-exp]   = {nstr(Theta_i_qexp, 20)}")

rel_err = fabs(Theta_i_qexp - Theta_i_formula) / fabs(Theta_i_formula)
log_err = float(log10(rel_err)) if rel_err > 0 else -999
print(f"  Relative error = {nstr(rel_err, 5)}  (log10 ≈ {log_err:.1f})")

check("rel_err(Theta_Leech q-exp vs formula) < 1e-50",
      rel_err < mpf('1e-50'),
      f"rel_err = {nstr(rel_err, 5)}")
check("Theta_Leech(i) > 0", Theta_i_formula > 0)
check("Theta_Leech(i) < 5", Theta_i_formula < 5)   # sanity: should be ~1.8

# Verify the formula components:
# (j(i) - 720) * Delta(i) = 1008 * eta(i)^24
formula_direct = mpf(1008) * Delta_i_cf
rel_err_direct = fabs(formula_direct - Theta_i_formula) / fabs(Theta_i_formula)
check("1008 * Delta(i) = (j(i)-720) * Delta(i)  [algebraic, < 1e-55]",
      rel_err_direct < mpf('1e-55'),
      f"rel_err = {nstr(rel_err_direct, 5)}")

# Verify 1008 = j(i) - 720
check("j(i) - 720 = 1728 - 720 = 1008", int(j_i) - 720 == 1008)
check("1008 = 2^4 * 3^2 * 7", 1008 == 2**4 * 3**2 * 7)

# Additional cross-checks
print("\n=== Additional cross-checks ===")

# The 24 coincidences
check("Leech dim = Niemeier count = 24", LEECH_DIM == 24)
check("Leech dim = eta exponent = 24", LEECH_DIM == ETA_EXPONENT)
check("Leech dim = B_G2/2 = 24", LEECH_DIM == int(B_G2) // 2)
check("Leech dim = 3*rank(E8) = 24", LEECH_DIM == 3 * RANK_E8)
check("Leech dim = Golay code length = 24", LEECH_DIM == GOLAY_LEN)

# J_short connections
check("Golay dim = J_short = 12", GOLAY_DIM == int(J_short))
check("j(i) = J_short^3 = 1728", int(j_i) == 1728)
check("12^3 = 1728", 12**3 == 1728)

# 720 = 3 * 240 = 3 * |Phi_E8|
check("720 = 3 * |Phi_E8| = 3 * 240", 720 == 3 * ROOTS_E8)
check("720 = [q^1] E4^3 (rootlessness correction)", e43_coeffs[1] == 720)

# Modular weight of Theta_Leech = 12 = Leech_dim / 2
check("Modular weight of Theta_Leech = 24/2 = 12", LEECH_DIM // 2 == 12)
check("Modular weight = J_short", LEECH_DIM // 2 == int(J_short))

# Theta_{Leech}(i) numerical sanity
theta_val = float(Theta_i_formula)
print(f"  Theta_Leech(i) ≈ {theta_val:.6f}  (expected ~1.8)")
check("1.5 < Theta_Leech(i) < 2.5  (sanity)", 1.5 < theta_val < 2.5)

# Co_0 order = 2 * Co_1 order
CO0_ORDER = 2 * CO1_ORDER
check("|Co_0| = 2^22 * 3^9 * 5^4 * 7^2 * 11 * 13 * 23",
      CO0_ORDER == 2**22 * 3**9 * 5**4 * 7**2 * 11 * 13 * 23)

# ── Summary ───────────────────────────────────────────────────────────────────
print("\n" + "="*65)
print("SUMMARY — Addendum P169 verification (The Leech Lattice)")
print("="*65)
print(f"  196560 = 2^4 * 3^3 * 5 * 7 * 13  ✓")
print(f"  196560 / 240 = {MIN_VECTORS // ROOTS_E8} = 3*273 = 9*7*13  ✓")
print(f"  Golay code dim = 12 = J_short  ✓")
print(f"  Leech dim = 24 = B_G2/2 = eta_exponent = 3*rank(E8)  ✓")
print(f"  [q^1] Theta_Leech = {leech_theta_coeffs[1]} (rootless)  ✓")
print(f"  [q^2] Theta_Leech = {leech_theta_coeffs[2]} (minimal vectors)  ✓")
print(f"  [q^3] Theta_Leech = {leech_theta_coeffs[3]}  ✓")
print(f"  [q^4] Theta_Leech = {leech_theta_coeffs[4]}  ✓")
print(f"  |Co_1| = {CO1_ORDER}  ✓")
print(f"  |Co_1| = 2^21*3^9*5^4*7^2*11*13*23  ✓")
print(f"  196560 * 12 = {MIN_VECTORS * 12}  ✓")
print(f"  Theta_Leech(i) [formula] = {nstr(Theta_i_formula, 15)}")
print(f"  Theta_Leech(i) [q-exp]   = {nstr(Theta_i_qexp, 15)}")
print(f"  rel_err = {nstr(rel_err, 5)}  (< 1e-50)  ✓")
print(f"  1008 * eta(i)^24 matches to 50 sig figs  ✓")
print(f"  j(i) - 720 = {int(j_i) - 720} = 1008  ✓")

print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
sys.exit(0 if FAIL == 0 else 1)