"""
verify_P159.py — Addendum P159: E₄(i)/η(i)^8 = J_short (closure of P158-OP1)

The j-function identity j(τ) = E₄(τ)³ / η(τ)²⁴  combined with
P158's result j(i) = J_short³ = 12³ directly forces E₄(i)/η(i)^8 = 12.
This script verifies the algebraic chain numerically at 50 d.p.

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

import sys
import mpmath
mpmath.mp.dps = 50

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}")

# Constants
J_short = 12
j_i = 1728  # = J_short^3 (P158)

# Compute E₄(i) and η(i) numerically at tau = i
tau = mpmath.mpc(0, 1)
q = mpmath.exp(2 * mpmath.pi * mpmath.j * tau)

# Dedekind eta: η(τ) = q^{1/24} ∏(1-q^n)
def eta(tau, terms=500):
    q = mpmath.exp(2 * mpmath.pi * mpmath.j * tau)
    result = q ** (mpmath.mpf(1)/24)
    for n in range(1, terms):
        result *= (1 - q**n)
    return result

# Eisenstein E₄: E₄(τ) = 1 + 240 Σ σ₃(n) q^n
def sigma3(n):
    return sum(d**3 for d in range(1, n+1) if n % d == 0)

def E4(tau, terms=200):
    q = mpmath.exp(2 * mpmath.pi * mpmath.j * tau)
    result = mpmath.mpf(1)
    for n in range(1, terms):
        result += 240 * sigma3(n) * q**n
    return result

eta_i = eta(tau)
E4_i = E4(tau)

# Check j = E4^3 / eta^24
j_computed = E4_i**3 / eta_i**24
print(f"j(i) computed = {float(mpmath.re(j_computed)):.6f}  (expect 1728)")

# Check x = E4/eta^8
x = E4_i / eta_i**8
print(f"E₄(i)/η(i)^8 = {float(mpmath.re(x)):.10f}  (expect 12)")
print(f"J_short = {J_short}")
print(f"j(i)^(1/3) = {j_i**(1/3):.10f}  (expect 12)")

# Verify chain
check(1, "j(i) = 1728", abs(mpmath.re(j_computed) - 1728) < 1e-6)
check(2, "E4/eta^8 = 12", abs(mpmath.re(x) - 12) < 1e-6)
check(3, "x = J_short", abs(mpmath.re(x) - J_short) < 1e-6)
print()
print("Chain:")
print(f"  j(τ) = E₄(τ)³/η(τ)²⁴  [definition]")
print(f"  j(i) = {j_i} = J_short³ = 12³  [P158]")
print(f"  E₄(i)/η(i)^8 = j(i)^{{1/3}} = J_short = {J_short}  [P159, QED]")
print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
sys.exit(0 if FAIL == 0 else 1)