#!/usr/bin/env python3
"""
verify_P230.py — Verifier for Addendum 230: Commutator Output Ratio -48*pi/(pi-1)
Checks a_2 = 3*pi^2*(1+8*pi)/16 exactly and fully factors the ratio
c_lin/c_const = -48*pi/(pi-1) = -4 * 12 * pi/(pi-1).

All arithmetic uses mpmath at dps=60 for exact verification.
"""

from mpmath import mp, mpf, pi as PI, fabs, nstr

mp.dps = 60

# ─── assertion harness ────────────────────────────────────────────────────────
PASS = 0; FAIL = 0
def check(name: str, condition: bool) -> None:
    global PASS, FAIL
    ok = bool(condition)
    if ok:
        PASS += 1
    else:
        FAIL += 1
    n = PASS + FAIL
    print(f"  [{'PASS' if ok else 'FAIL'}] {n:>2}. {name}")

OMEGA_0 = PI**3 / 4     # = a_3 (A225)
OMEGA   = 4*PI**3 + PI**2 + PI

# ─── Section 1: Exact a_2 ────────────────────────────────────────────────────
print("\nS1  Exact a_2 = 3*pi^2*(1+8*pi)/16")

# Derivation: [coeff of x^2 in rho] = 3*pi^2
# From Gegenbauer inversion: 3*pi^2 = -96*a_3 + 16*a_2
# => a_2 = (3*pi^2 + 96*a_3) / 16 = (3*pi^2 + 24*pi^3) / 16

a3 = OMEGA_0
a2_formula = (3*PI**2 + 24*PI**3) / 16
a2_closed  = 3 * PI**2 * (1 + 8*PI) / 16

check("P001  a_2 = (3*pi^2 + 24*pi^3)/16  (from x^2 coeff inversion)",
      fabs(a2_formula - a2_closed) < mpf('1e-55'))

check("P002  a_2 = 3*pi^2*(1+8*pi)/16  (factored form)",
      fabs(a2_closed - (3*PI**2*(1+8*PI)/16)) < mpf('1e-55'))

check("P003  a_2 > 0  (positive)",
      a2_closed > 0)

# Numerical value approx 48.36
check("P004  a_2 in [48, 49]  (sanity bound)",
      48 < float(a2_closed) < 49)

# Cross-check: x^2 coefficient of rho reconstructed
x2_coeff_rho = 3 * PI**2
reconstructed_x2 = -96*a3 + 16*a2_closed
check("P005  -96*a_3 + 16*a_2 = 3*pi^2  (Gegenbauer roundtrip)",
      fabs(reconstructed_x2 - x2_coeff_rho) < mpf('1e-55'))

# Bulk vs boundary decomposition of a_2:
# a_2 = 3*pi^2/16 * (1 + 8*pi) = (3*pi^2/16) + (3*pi^2/16)*8*pi
a2_boundary_part = 3*PI**2 / 16
a2_bulk_correction = 3*PI**2 * 8*PI / 16

check("P006  a_2 = boundary_part + bulk_correction where boundary = 3*pi^2/16",
      fabs(a2_boundary_part + a2_bulk_correction - a2_closed) < mpf('1e-55'))

check("P007  bulk_correction / boundary_part = 8*pi",
      fabs((a2_bulk_correction / a2_boundary_part) - 8*PI) < mpf('1e-55'))

# ─── Section 2: Commutator Gegenbauer coefficients ────────────────────────────
print("\nS2  C_0, C_1 in Gegenbauer basis")

# From A226 Theorem 3:
# [D^2,Delta]rho = (-24576*a_2 - 276480*Omega_0)*U_0 + (-221184*Omega_0)*U_1
C1 = -221184 * OMEGA_0
C0 = -24576 * a2_closed - 276480 * OMEGA_0

check("P008  C_1 = -221184 * Omega_0 = -221184 * pi^3/4",
      fabs(C1 - (-221184 * PI**3 / 4)) < mpf('1e-55'))

check("P009  C_1 = -55296 * pi^3",
      fabs(C1 - (-55296 * PI**3)) < mpf('1e-55'))

# C_0 exact computation:
# = -24576*(3*pi^2*(1+8*pi)/16) - 276480*(pi^3/4)
# = -4608*pi^2*(1+8*pi) - 69120*pi^3
# = -4608*pi^2 - 36864*pi^3 - 69120*pi^3
# = -4608*pi^2 - 105984*pi^3
C0_exact = -mpf('4608')*PI**2 - mpf('105984')*PI**3

check("P010  C_0 = -4608*pi^2 - 105984*pi^3  (Gegenbauer scalar coeff)",
      fabs(C0 - C0_exact) < mpf('1e-50'))

# ─── Section 3: Monomial output ───────────────────────────────────────────────
print("\nS3  Monomial form c_const + c_lin*x")

# c_const = C_0 - 2*C_1
# c_lin   = 4*C_1
c_const = C0 - 2*C1
c_lin   = 4*C1

# Expected: c_const = 4608*pi^2*(pi-1), c_lin = -221184*pi^3
c_const_expected = mpf('4608') * PI**2 * (PI - 1)
c_lin_expected   = -mpf('221184') * PI**3

check("P011  c_const = C_0 - 2*C_1 = 4608*pi^2*(pi-1)",
      fabs(c_const - c_const_expected) < mpf('1e-50'))

check("P012  c_lin = 4*C_1 = -221184*pi^3",
      fabs(c_lin - c_lin_expected) < mpf('1e-50'))

# Algebraic derivation of c_const:
# C_0 - 2*C_1 = (-4608*pi^2 - 105984*pi^3) - 2*(-55296*pi^3)
#             = -4608*pi^2 - 105984*pi^3 + 110592*pi^3
#             = -4608*pi^2 + 4608*pi^3
#             = 4608*pi^2*(pi - 1)
step1 = -mpf('4608')*PI**2 - mpf('105984')*PI**3
step2 = step1 + 2*mpf('55296')*PI**3   # subtract 2*C1
step3 = -mpf('4608')*PI**2 + mpf('4608')*PI**3

check("P013  Derivation step: -105984 + 110592 = +4608  (bulk cancellation)",
      fabs(mpf('-105984') + mpf('110592') - 4608) < mpf('1e-55'))

check("P014  c_const = -4608*pi^2 + 4608*pi^3 = 4608*pi^2*(pi-1)  (algebraic)",
      fabs(step3 - c_const_expected) < mpf('1e-55'))

# Sign check
check("P015  c_const > 0  (since pi > 1, pi-1 > 0)",
      c_const > 0)

check("P016  c_lin < 0  (negative linear component)",
      c_lin < 0)

# ─── Section 4: The ratio -48*pi/(pi-1) ──────────────────────────────────────
print("\nS4  Ratio c_lin/c_const = -48*pi/(pi-1)")

ratio = c_lin / c_const
ratio_formula = -48 * PI / (PI - 1)

check("P017  c_lin/c_const = -48*pi/(pi-1)  (exact closed form)",
      fabs(ratio - ratio_formula) < mpf('1e-50'))

check("P018  Ratio is negative (linear term dominates sign)",
      ratio < 0)

check("P019  |ratio| approx 70.4  (sanity bound: 69 < |r| < 72)",
      69 < float(fabs(ratio)) < 72)

# ─── Section 5: Factorisation of 48 ─────────────────────────────────────────
print("\nS5  48 = ell_1 * g_{31}")

ell_1 = mpf('4')    # leading coeff of U_1 = 4x-2 in Hopf x-coordinate
g_31  = mpf('12')   # |mu_3 - mu_1| = |-15 - (-3)| = 12

check("P020  ell_1 = 4  (leading coefficient of U_1 = 4x-2)",
      fabs(ell_1 - 4) < mpf('1e-55'))

check("P021  g_{31} = |mu_3 - mu_1| = |-15 - (-3)| = 12",
      fabs(g_31 - 12) < mpf('1e-55'))

check("P022  48 = ell_1 * g_{31} = 4 * 12  (integer factorisation)",
      fabs(ell_1 * g_31 - 48) < mpf('1e-55'))

check("P023  -221184 / 4608 = -48  (exact integer quotient)",
      fabs(mpf('-221184') / mpf('4608') - (-48)) < mpf('1e-55'))

# Cross-check: 221184 = ell_1 * g_31 * 4608
check("P024  221184 = 48 * 4608  (spectral integers * K_4)",
      fabs(48 * mpf('4608') - 221184) < mpf('1e-55'))

check("P025  4608 = K_{d=4}  (prefactor from A227)",
      fabs(mpf('4608') - 128*1*2*3*6) < mpf('1e-55'))

# ─── Section 6: pi/(pi-1) — layer arithmetic ─────────────────────────────────
print("\nS6  pi/(pi-1) — layer fraction analysis")

pi_ratio = PI / (PI - 1)

check("P026  pi/(pi-1) is the transcendental factor in the ratio",
      fabs(ratio_formula - (-48 * pi_ratio)) < mpf('1e-55'))

check("P027  pi > 1  (so pi-1 > 0 and pi/(pi-1) > 0)",
      PI > 1)

# Layer fractions
f_edge = PI / OMEGA
f_bnd  = PI**2 / OMEGA
f_bulk = 4*PI**3 / OMEGA

check("P028  f_edge + f_bnd + f_bulk = 1.0  (partition of unity)",
      fabs(f_edge + f_bnd + f_bulk - 1) < mpf('1e-55'))

# pi/(pi-1) != any simple combination of layer fractions (open sub-question)
# Check that pi/(pi-1) != f_bulk/f_bnd etc.
check("P029  pi/(pi-1) != f_bulk/f_edge  (not a simple inter-layer ratio)",
      fabs(pi_ratio - f_bulk/f_edge) > mpf('1'))

check("P030  pi/(pi-1) in (1.0, 2.0)  (numerical range)",
      1 < float(pi_ratio) < 2)

# The (pi-1) factor traces to bulk-boundary cancellation in a_2:
# -24576 * a_2 = -24576 * 3*pi^2*(1+8*pi)/16 = -4608*pi^2 - 36864*pi^3
# -276480 * Omega_0 = -69120*pi^3
# Sum = -4608*pi^2 - (36864+69120)*pi^3 = -4608*pi^2 - 105984*pi^3 = C_0
# Then C_0 - 2*C_1 adds back 2*55296*pi^3 = 110592*pi^3
# Net pi^3 coeff: -105984 + 110592 = +4608
# Net pi^2 coeff: -4608
# = 4608*(pi^3 - pi^2) = 4608*pi^2*(pi-1)

check("P031  -105984 + 110592 = 4608  (bulk-boundary residue at n=2)",
      fabs(mpf('-105984') + mpf('110592') - 4608) < mpf('1e-55'))

check("P032  4608*(pi^3 - pi^2) = 4608*pi^2*(pi-1)  (factored form)",
      fabs(mpf('4608')*(PI**3 - PI**2) - c_const_expected) < mpf('1e-55'))

# ─── Section 7: Numerical consistency ────────────────────────────────────────
print("\nS7  Numerical consistency (mp.dps=60)")

check("P033  c_lin = -221184 * pi^3  (high-precision)",
      fabs(c_lin + 221184*PI**3) < mpf('1e-50'))

check("P034  c_const = 4608 * pi^2 * (pi-1)  (high-precision)",
      fabs(c_const - 4608*PI**2*(PI-1)) < mpf('1e-50'))

check("P035  ratio = c_lin/c_const = -48*pi/(pi-1)  (high-precision)",
      fabs(ratio + 48*PI/(PI-1)) < mpf('1e-50'))

# Additional: verify 4608 * pi^2 * (pi-1) is positive and order ~300
check("P036  c_const in [9e4, 1.1e5]  (numerical sanity: 4608*pi^2*(pi-1) ~97400)",
      9e4 < float(c_const) < 1.1e5)

# 221184*pi^3 approx 221184 * 31.006 ~ 6.86e6
check("P037  |c_lin| in [5e6, 8e6]  (numerical sanity)",
      5e6 < float(fabs(c_lin)) < 8e6)

# ─── Section 8: Integer crosschecks ──────────────────────────────────────────
print("\nS8  Integer structure")

check("P038  221184 = 18432 * 12  (equal-weight * spectral gap)",
      fabs(mpf('221184') - 18432*12) < mpf('1e-55'))

check("P039  276480 = 18432 * 15  (entry [0,3] in commutator * |mu_3|)",
      fabs(mpf('276480') - 18432*15) < mpf('1e-55'))

check("P040  24576 = 3072 * 8  (entry [0,2] in D^2_geg * |mu_2 - mu_0|)",
      fabs(mpf('24576') - 3072*8) < mpf('1e-55'))

# Summary
print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
if FAIL:
    raise SystemExit(1)