"""
verify_P157.py — Addendum P157: TOE constants and modular forms at CM points

Systematic search: does the TOE constant structure appear in modular forms
at CM points?

Four computations:
  1. j(τ_D) at CM points for fundamental discriminants D ∈ {3..163}
  2. Eta and Eisenstein series at τ = i
  3. Ramanujan near-integer scan: e^{π√D} vs α^k for D ∈ [1,300]
  4. Depth of the 1728 connection

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

import math
import sys
import mpmath

mpmath.mp.dps = 60   # 60 decimal places

PASS = FAIL = 0
_N = 0

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

# ── TOE constants ──────────────────────────────────────────────────────────────
ALPHA_INV  = 4*math.pi**3 + math.pi**2 + math.pi   # ≈ 137.036
BREATH     = math.pi * ALPHA_INV                     # ≈ 430.512
PHI        = (1 + 5**0.5) / 2                        # ≈ 1.618
SIN2_W     = 3 / (8 * PHI)                           # ≈ 0.23176
B_G2       = 48      # Killing form G₂ on H_s (P151)
B_F4       = 36      # Killing form F₄ on H_s (P153)
J_SHORT    = 12      # short-root sector
J_LONG     = 36      # long-root sector
PRODUCT    = 432     # J_short × J_long
CUBE       = 1728    # B_F4 × B_G2 = 12³ = j(i)
TOL        = 0.005   # 0.5% tolerance

TOE_TARGETS = {
    "ALPHA_INV": ALPHA_INV, "BREATH": BREATH, "PHI": PHI,
    "SIN2_W"   : SIN2_W,   "B_G2"  : B_G2,   "B_F4": B_F4,
    "J_SHORT"  : J_SHORT,  "J_LONG": J_LONG,
    "PRODUCT"  : PRODUCT,  "CUBE"  : CUBE,
}

def near(val, target, tol=TOL):
    if target == 0: return abs(val) < tol
    return abs(float(val)/float(target) - 1) < tol

print("=" * 72)
print("P157 VERIFICATION: TOE CONSTANTS vs MODULAR FORMS AT CM POINTS")
print("=" * 72)
print(f"ALPHA_INV  = {ALPHA_INV:.10f}")
print(f"BREATH     = {BREATH:.10f}")
print(f"PRODUCT    = {PRODUCT}  (= 12×36 = J_short×J_long)")
print(f"CUBE       = {CUBE}  (= 12³ = B_F4×B_G2 = j(i))")
print(f"Tolerance  = {TOL*100:.1f}%")
print()

# ── Utility: σ_k(n) ───────────────────────────────────────────────────────────
def sigma(k, n):
    return sum(d**k for d in range(1, n+1) if n % d == 0)

# ── Utility: Eisenstein E4, E6 via q-expansion ────────────────────────────────
def eisenstein(tau, terms=350):
    q = mpmath.exp(2 * mpmath.pi * mpmath.j * tau)
    E4 = mpmath.mpf(1)
    E6 = mpmath.mpf(1)
    qn = q
    for n in range(1, terms+1):
        E4 += 240 * sigma(3, n) * qn
        E6 -= 504 * sigma(5, n) * qn
        qn *= q
        if abs(qn) < mpmath.mpf(10)**(-55):
            break
    return E4, E6

# ── Utility: η(τ) via q-product ───────────────────────────────────────────────
def eta_func(tau, terms=2000):
    q = mpmath.exp(2 * mpmath.pi * mpmath.j * tau)
    prod = q**(mpmath.mpf(1)/24)
    for n in range(1, terms+1):
        prod *= (1 - q**n)
        if abs(q**n) < mpmath.mpf(10)**(-55):
            break
    return prod

# ══════════════════════════════════════════════════════════════════════════════
# COMPUTATION 1 — j(τ_D) at CM points
# ══════════════════════════════════════════════════════════════════════════════
print("─" * 72)
print("COMPUTATION 1: j(τ_D) AT CM POINTS")
print("─" * 72)
print()

# Exact integer j-values (classical CM theory)
cm_j_exact = {
    3:   0,           7:   -3375,       4:   1728,
    8:   8000,        11:  -32768,      12:  54000,
    16:  287496,      19:  -884736,     27:  -12288000,
    28:  16581375,    43:  -884736000,  67:  -147197952000,
    163: -262537412640768000,
}

# CM points τ_D
def cm_tau(D):
    sqrt_D = mpmath.sqrt(D)
    if D in [3, 7, 11, 19, 27, 43, 67, 163]:
        return mpmath.mpc(mpmath.mpf(1)/2, sqrt_D/2)
    else:  # D=4→i, D=8→i√2, D=12→i√3, D=16→2i, D=28→i√7
        return mpmath.mpc(0, sqrt_D/2 if D in [4, 8, 12] else sqrt_D/mpmath.sqrt(4))

cm_tau_map = {
    3:  mpmath.mpc(mpmath.mpf(1)/2, mpmath.sqrt(3)/2),
    4:  mpmath.mpc(0, 1),
    7:  mpmath.mpc(mpmath.mpf(1)/2, mpmath.sqrt(7)/2),
    8:  mpmath.mpc(0, mpmath.sqrt(2)),
    11: mpmath.mpc(mpmath.mpf(1)/2, mpmath.sqrt(11)/2),
    12: mpmath.mpc(0, mpmath.sqrt(3)),
    16: mpmath.mpc(0, 2),
    19: mpmath.mpc(mpmath.mpf(1)/2, mpmath.sqrt(19)/2),
    27: mpmath.mpc(mpmath.mpf(1)/2, mpmath.sqrt(27)/2),
    28: mpmath.mpc(0, mpmath.sqrt(7)),
    43: mpmath.mpc(mpmath.mpf(1)/2, mpmath.sqrt(43)/2),
    67: mpmath.mpc(mpmath.mpf(1)/2, mpmath.sqrt(67)/2),
    163:mpmath.mpc(mpmath.mpf(1)/2, mpmath.sqrt(163)/2),
}

print(f"{'D':>4}  {'j(τ_D) exact':>25}  {'q-exp verified':>5}  {'TOE hits on |j| or |j|^(1/3)'}")
print("-" * 85)

cm_flags = {}
for D in sorted(cm_j_exact):
    j_ex = cm_j_exact[D]
    tau_mp = cm_tau_map[D]
    E4, E6 = eisenstein(tau_mp, terms=200 if D < 20 else 80)
    Delta = (E4**3 - E6**2)/1728
    j_qexp = float(mpmath.re(E4**3 / Delta)) if abs(Delta) > 1e-300 else 0.0
    ok = abs(j_qexp - j_ex) < 2.0

    hits = []
    abs_j = abs(j_ex)
    if abs_j > 0:
        cbrt_j = abs_j**(1/3)
        for name, t in TOE_TARGETS.items():
            if t > 0 and near(abs_j, t): hits.append(f"|j|={name}")
            if t > 0 and near(cbrt_j, t): hits.append(f"|j|^(1/3)={name}")
            if t > 0 and near(abs_j/4, t): hits.append(f"|j|/4={name}")
    cm_flags[D] = hits
    hit_str = "; ".join(hits) if hits else "—"
    check(f"D={D}: j(τ_D) exact {j_ex} verified by q-expansion; TOE hits: {hit_str}", ok)

print()
cbrt_27 = abs(cm_j_exact[27])**(1/3)
cbrt_43 = abs(cm_j_exact[43])**(1/3)
cbrt_163 = abs(cm_j_exact[163])**(1/3)
print(f"  |j(τ_27)|^(1/3) = {cbrt_27:.6f}  (not a TOE constant; ≠ 12k for integer k)")
print(f"  |j(τ_43)|^(1/3) = {cbrt_43:.6f}  (= 960 exactly; 960/ALPHA_INV = {cbrt_43/ALPHA_INV:.4f})")
print(f"  |j(τ_163)|^(1/3) = {cbrt_163:.6f}  (= 640320)")
print(f"  Only D=4 yields j(i)=1728=12³ with 12=J_short  ← the P153 origin")

# ══════════════════════════════════════════════════════════════════════════════
# COMPUTATION 2 — Eta and Eisenstein at τ = i
# ══════════════════════════════════════════════════════════════════════════════
print()
print("─" * 72)
print("COMPUTATION 2: ETA AND EISENSTEIN SERIES AT τ = i")
print("─" * 72)
print()

tau_i = mpmath.mpc(0, 1)
gamma_14 = mpmath.gamma(mpmath.mpf(1)/4)

# η(i)
eta_i = eta_func(tau_i)
eta_i_known = gamma_14 / (2 * mpmath.pi**(mpmath.mpf(3)/4))
print(f"η(i) via q-product     = {float(abs(eta_i)):.20f}")
print(f"η(i) = Γ(1/4)/(2π^3/4) = {float(abs(eta_i_known)):.20f}")
check("η(i) q-product matches Γ(1/4)/(2π^3/4) (60 dps)",
      mpmath.almosteq(abs(eta_i), abs(eta_i_known), 1e-30))
print()

# E4(i), E6(i)
E4_i, E6_i = eisenstein(tau_i, terms=350)
Delta_i_alt = (E4_i**3 - E6_i**2)/1728
j_i_num = float(mpmath.re(E4_i**3 / Delta_i_alt))
E4_i_re = float(mpmath.re(E4_i))
E6_i_re = float(mpmath.re(E6_i))

print(f"E4(i)                  = {E4_i_re:.20f}")
print(f"E6(i)                  = {E6_i_re:.4e}  (→ 0 by symmetry of CM(i))")
print(f"j(i) = E4³/Δ           = {j_i_num:.6f}  (exact: 1728)")
print()

# Known exact: E4(i) = 3Γ(1/4)^8 / (2π)^6
E4_known = 3 * gamma_14**8 / (2*mpmath.pi)**6
print(f"E4(i) known formula 3Γ(1/4)^8/(2π)^6 = {float(mpmath.re(E4_known)):.20f}")
check("E4(i) q-expansion matches known formula 3Γ(1/4)^8/(2π)^6",
      mpmath.almosteq(E4_i, E4_known, 1e-20))
print()

# KEY FINDING: E4(i) / η(i)^8
eta_i_8 = eta_i**8
ratio_E4_eta8 = float(mpmath.re(E4_i / eta_i_8))
print(f"*** E4(i) / η(i)^8     = {ratio_E4_eta8:.10f}  ***")
print(f"    This equals J_SHORT = 12 exactly.")
print(f"    Proof: η(i)^8 = Γ(1/4)^8/(2^8 π^6);  E4(i) = 3Γ(1/4)^8/(2π)^6")
print(f"           E4(i)/η(i)^8 = 3/(2π)^6 × (2π^(3/4))^8 = 3×2^2 = 12")
print()

# η(i)/η(2i) Hauptmodul check
eta_2i = eta_func(mpmath.mpc(0, 2))
t2_i = float(mpmath.re((eta_i/eta_2i)**24))
print(f"(η(i)/η(2i))^24        = {t2_i:.6f}  (= 2^9 = 512, not a TOE target)")
print()

# Ratios of E4(i) to TOE constants
print("E4(i) ratios (no near-match within 0.5% found):")
print(f"  E4(i)/ALPHA_INV      = {E4_i_re/ALPHA_INV:.8f}")
print(f"  E4(i)/BREATH         = {E4_i_re/BREATH:.8f}")
print(f"  E4(i)/PRODUCT        = {E4_i_re/PRODUCT:.8f}")
print(f"  E4(i)^(1/4)          = {E4_i_re**(1/4):.8f}")
E8_i_re  = E4_i_re**2
E12_i_re = (441*E4_i_re**3 + 250*E6_i_re**2)/691
print(f"  E8(i)=E4(i)²         = {E8_i_re:.8f}")
print(f"  E12(i)               = {E12_i_re:.8f}")
print("  (No Eisenstein ratio to π^k or ALPHA_INV within 0.5%)")

# ══════════════════════════════════════════════════════════════════════════════
# COMPUTATION 3 — Ramanujan near-integer scan
# ══════════════════════════════════════════════════════════════════════════════
print()
print("─" * 72)
print("COMPUTATION 3: RAMANUJAN NEAR-INTEGER SCAN  e^{π√D} vs α^{-k}")
print("─" * 72)
print()
print("Threshold: |ratio − nearest_int| < 0.01  (1%)")
print("PRECISION NOTE: for D ≥ 43 with k=±3, e^{π√D}/α^3 exceeds 10^{14},")
print("  losing fractional precision in float64 → spurious error≈0. Meaningful")
print("  regime is D ≤ 40 (float64 reliable to 15 sig figs).")
print()

mpmath.mp.dps = 30
print("D ≤ 40 hits (float64-reliable regime):")
small_hits = []
for D in range(1, 41):
    val = float(mpmath.exp(mpmath.pi * mpmath.sqrt(D)))
    for k in [1, 2, 3, -1, -2, -3]:
        denom = ALPHA_INV**k
        ratio = val / denom
        nearest_int = round(ratio)
        if nearest_int > 0:
            err = abs(ratio - nearest_int)
            if err < 0.01:
                small_hits.append((D, k, nearest_int, err))
                print(f"  D={D:2d}: e^{{π√{D}}} ≈ {nearest_int} × α^{{{-k}}},"
                      f"  error={err:.6f},  e^{{π√{D}}}={val:.2f}")

expected_random = 2 * 6 * 40 * 0.01   # 2 sides × 6 powers × 40 D × 1% threshold
print(f"\nHits in [1,40]: {len(small_hits)}  (expected from randomness: ~{expected_random:.1f})")
print("None of these integers has structural significance in the TOE.")
print()
print("D = 43 (Heegner): e^{π√43} ≈ 884736744 (= |j(τ_43)| + 744)")
val43 = float(mpmath.exp(mpmath.pi * mpmath.sqrt(43)))
print(f"  e^{{π√43}} = {val43:.4f}")
print(f"  |j(τ_43)| + 744 = {abs(cm_j_exact[43]) + 744}")
print(f"  e^{{π√43}} / ALPHA_INV^3 = {val43 / ALPHA_INV**3:.6f}  (not near integer)")
mpmath.mp.dps = 60

# ══════════════════════════════════════════════════════════════════════════════
# COMPUTATION 4 — Depth of the 1728 connection
# ══════════════════════════════════════════════════════════════════════════════
print()
print("─" * 72)
print("COMPUTATION 4: DEPTH OF THE 1728 CONNECTION")
print("─" * 72)
print()

G1_abs = PRODUCT - BREATH   # 432 - π·α⁻¹

print("Exact integer arithmetic:")
check("1728 = 12³", 1728 == 12**3)
check("1728 = B_F4 × B_G2 = 36×48", 1728 == B_F4 * B_G2)
print(f"  1728 = j(i)                  classical CM theorem")
check("1728 / 4 = 432 = PRODUCT", 1728//4 == PRODUCT)
check("1728^(1/3) = 12 = J_SHORT", round(1728**(1/3)) == J_SHORT)
print(f"  12 appears: J_SHORT, ∛j(i), E4(i)/η(i)^8 — three independent routes")
print()
print("BREATH vs j(i)/4:")
print(f"  j(i)/4 = 1728/4 = 432.000000...  (exact integer = PRODUCT)")
print(f"  BREATH = π·α⁻¹ = {BREATH:.12f}")
print(f"  j(i)/4 − BREATH = G1_abs    = {G1_abs:.12f}  (P156, the known unknown)")
print(f"  relative gap    = G1_rel    = {G1_abs/BREATH:.12f}  ({G1_abs/BREATH*100:.4f}%)")
print(f"  j(i)/4 is PRODUCT, not BREATH; the G1 gap separates them.")
print()
print("Structural summary at τ=i:")
print(f"  j(i)          = 1728 = 12³ = B_F4×B_G2")
print(f"  j(i)^(1/3)    = 12   = J_SHORT")
print(f"  j(i)/4        = 432  = PRODUCT = J_SHORT × J_LONG")
print(f"  E4(i)/η(i)^8  = 12   = J_SHORT  (provable identity)")
print(f"  E6(i)         = 0    (CM symmetry at τ=i, order 4)")
print(f"  BREATH        = {BREATH:.6f} ≠ 432  (gap = G1_abs)")

# ══════════════════════════════════════════════════════════════════════════════
# SUMMARY AND VERDICT
# ══════════════════════════════════════════════════════════════════════════════
print()
print("─" * 72)
print("SUMMARY OF NEAR-MATCHES (TOL = 0.5%)")
print("─" * 72)
print()
print("CM TABLE (D=3..163):")
for D, hits in sorted(cm_flags.items()):
    if hits:
        print(f"  D={D}: {'; '.join(hits)}")
    else:
        print(f"  D={D}: no match")

print()
print("EISENSTEIN at τ=i:")
print("  E4(i)/η(i)^8 = 12.0000000000 = J_SHORT  (exact classical identity)")
print("  No other Eisenstein ratio within 0.5% of any TOE target")

print()
print("RAMANUJAN D≤40:")
if small_hits:
    for D, k, ni, err in small_hits:
        print(f"  D={D}: error={err:.4f} (within statistical expectation)")
else:
    print("  No hits")

print()
print("─" * 72)
print("VERDICT")
print("─" * 72)
print("""
1. j(i)=1728=B_{F₄}×B_{G₂}: one real coincidence, both sides = 1728 exactly.
   D=4 is the only discriminant where any CM j-value touches a TOE constant.

2. E4(i)/η(i)^8 = 12 exactly. This is a classical provable identity; the
   coincidence with J_SHORT=12 is real but 12 is already determined by
   j(i)=12³. One underlying fact (j(i)=1728), three appearances of 12.

3. j(i)/4 = 432 = PRODUCT (exact integers). BREATH = π·α⁻¹ ≈ 430.512.
   The gap is G1_abs (P156) — the documented known unknown.

4. Ramanujan scan D≤40: hits consistent with statistical expectation.
   No hit approaches the structural necessity of the Heegner-163 result.
   D≥43 results are float64 precision artifacts, not mathematical truths.

CONCLUSION: j(i)=1728 is a genuine, isolated structural contact between CM
theory and TOE root-system geometry. It does not extend to other discriminants
or modular forms. The number 12 (=J_SHORT=∛j(i)=E4(i)/η(i)^8) is the common
thread, but its appearance in classical modular theory is independent of the
TOE derivation. Evidence is consistent with coincidence at D=4 — significant
but isolated, not systematic.
""")
print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
sys.exit(0 if FAIL == 0 else 1)
