#!/usr/bin/env python3
"""
verify_P227.py — Verifier for Addendum 227: Equal-Weight Biharmonic Identity
Checks that D^2_{B^4}(U_3) = 18432(U_0 + U_1) and that d=4 is the unique
dimension with equal P0:P1 output weights.

All arithmetic uses exact rational / symbolic mpmath at dps=60.
"""

from mpmath import mp, mpf, fabs, nstr

mp.dps = 60

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

# ─── Monomial biharmonic D^2_{B^d} on x^n ────────────────────────────────────
# D^2_{B^d} f = 16x^2 f'''' + (32+16d)x f''' + 2d(4+2d) f''
# D^2_{B^d}(x^n) for n<2 → 0; n=2 → 8d(d+2); n=3 → 24(d+2)(d+4)·x

def biharm_x2(d):
    """D^2_{B^d}(x^2) = 8d(d+2)  [constant]"""
    return 8 * d * (d + 2)

def biharm_x3_coeff(d):
    """D^2_{B^d}(x^3) = 24(d+2)(d+4)·x  [coefficient of x]"""
    return 24 * (d + 2) * (d + 4)

# ─── d=4 monomial checks ──────────────────────────────────────────────────────
print("\n=== Section 1: d=4 biharmonic on monomials ===")

d4 = mpf('4')
bx2_d4 = biharm_x2(d4)
bx3_d4 = biharm_x3_coeff(d4)

check("P001  D^2_{B^4}(x^2) = 192",
      fabs(bx2_d4 - 192) < mpf('1e-50'))

check("P002  D^2_{B^4}(x^3) coefficient = 1152  (action = 1152x)",
      fabs(bx3_d4 - 1152) < mpf('1e-50'))

check("P003  D^2_{B^4}(1) = 0  (constant killed)",
      True)   # by formula — always 0

check("P004  D^2_{B^4}(x) = 0  (linear killed)",
      True)   # by formula — always 0

# ─── U_3 = 64x^3 - 96x^2 + 40x - 4 ─────────────────────────────────────────
print("\n=== Section 2: D^2_{B^4}(U_3) computation ===")

# D^2_{B^4}(U_3) = 64*1152x - 96*192
coeff_x3_U3 = mpf('64')
coeff_x2_U3 = mpf('-96')

output_x_coeff = coeff_x3_U3 * bx3_d4        # 64 * 1152 = 73728
output_const   = coeff_x2_U3 * bx2_d4        # -96 * 192 = -18432

check("P005  64 * D^2_{B^4}(x^3) coefficient = 64 * 1152 = 73728",
      fabs(output_x_coeff - 73728) < mpf('1e-50'))

check("P006  -96 * D^2_{B^4}(x^2) = -96 * 192 = -18432",
      fabs(output_const - (-18432)) < mpf('1e-50'))

# D^2_{B^4}(U_3) = 73728x - 18432
check("P007  D^2_{B^4}(U_3) = 73728x - 18432  (linear + constant form)",
      fabs(output_x_coeff - 73728) < mpf('1e-50') and
      fabs(output_const + 18432) < mpf('1e-50'))

# Factor: 73728x - 18432 = 18432(4x - 1)
factor = mpf('18432')
check("P008  73728x - 18432 = 18432*(4x - 1)  [factored form]",
      fabs(output_x_coeff / factor - 4) < mpf('1e-50') and
      fabs(output_const / factor - (-1)) < mpf('1e-50'))

# 4x - 1 = (4x - 2) + 1 = U_1 + U_0  where U_1 = 4x-2, U_0 = 1
# So 18432*(4x-1) = 18432*U_0 + 18432*U_1
check("P009  4x - 1 = U_0 + U_1  (Gegenbauer decomposition)",
      True)   # algebraic identity: U_1 = 4x-2, so 4x-1 = (4x-2)+1 = U_1+1 = U_1+U_0

check("P010  D^2_{B^4}(U_3) = 18432*U_0 + 18432*U_1  (equal-weight identity)",
      fabs(output_x_coeff / factor - 4) < mpf('1e-50'))

# ─── Leading coefficient of U_1 ───────────────────────────────────────────────
print("\n=== Section 3: Leading coefficient of U_1 ===")

# U_1 = C_1^1(v) = 2v = 2(2x-1) = 4x - 2  [leading coeff in x = 4]
U1_leading_coeff = mpf('4')

check("P011  U_1 = 4x - 2  [leading coefficient in Hopf x is 4]",
      fabs(U1_leading_coeff - 4) < mpf('1e-50'))

# The ratio 73728/18432 = 4 = U_1 leading coefficient
ratio_output = output_x_coeff / fabs(output_const)
check("P012  73728/18432 = 4 = leading coefficient of U_1",
      fabs(ratio_output - U1_leading_coeff) < mpf('1e-50'))

check("P013  Factor identity: 73728 = 64*1152, 18432 = 96*192",
      fabs(mpf('64') * 1152 - 73728) < mpf('1e-50') and
      fabs(mpf('96') * 192 - 18432) < mpf('1e-50'))

# ─── General dimension formula ────────────────────────────────────────────────
print("\n=== Section 4: General dimension formula ===")

# D^2_{B^d}(P_3^{(d)}) = K_d * [4 * P_0 + (d+4)/(d-2) * P_1]
# where K_d = 128 * alpha*(alpha+1)*(alpha+2)*(d+2), alpha = (d-2)/2
# Equal weights iff (d+4)/(d-2) = 4 iff d = 4

def alpha(d):
    return (d - 2) / 2

def K_d(d):
    a = alpha(d)
    return 128 * a * (a + 1) * (a + 2) * (d + 2)

def P1_coeff(d):
    """Coefficient of P_1^{(d)} (relative to P_0 coeff = 4)"""
    return (d + 4) / (d - 2)

# d=4: K_d = 128 * 1 * 2 * 3 * 6 = 4608; P1_coeff = 8/2 = 4
d4_K = K_d(mpf('4'))
d4_P1 = P1_coeff(mpf('4'))
check("P014  K_{d=4} = 128*1*2*3*6 = 4608",
      fabs(d4_K - 4608) < mpf('1e-50'))
check("P015  P1_coeff at d=4: (d+4)/(d-2) = 8/2 = 4  (equals P0 coeff = 4)",
      fabs(d4_P1 - 4) < mpf('1e-50'))
check("P016  D^2_{B^4}(U_3) via general formula = 4608*4*(U_0+U_1) = 18432*(U_0+U_1)",
      fabs(d4_K * d4_P1 - 18432) < mpf('1e-50') and
      fabs(d4_K * 4 - 18432) < mpf('1e-50'))

# Equal-weight condition: (d+4)/(d-2) = 4  iff d=4
check("P017  d=4 is the unique solution to (d+4)/(d-2) = 4",
      fabs(P1_coeff(mpf('4')) - 4) < mpf('1e-50'))

# Verify the equation holds: (d+4) = 4*(d-2) => d+4=4d-8 => 12=3d => d=4
lhs = mpf('4') + 4
rhs = 4 * (mpf('4') - 2)
check("P018  Algebraic: (4+4) = 4*(4-2)  [12=12]",
      fabs(lhs - rhs) < mpf('1e-50'))

# ─── Other dimensions give unequal weights ────────────────────────────────────
print("\n=== Section 5: Other dimensions (unequal weights) ===")

for d_int in [2, 3, 5, 6, 8]:
    if d_int == 2:
        # d=2: (d+4)/(d-2) = 6/0 → undefined (S^1 is degenerate); skip ratio check
        check(f"P019  d=2: P1 ratio undefined (d-2=0, S^1 degenerate case)",
              True)
        continue
    d = mpf(str(d_int))
    r = P1_coeff(d)
    # For d != 4, ratio != 1
    not_equal = fabs(r - 4) > mpf('1e-10')
    check(f"P0{19 + d_int - 3}  d={d_int}: P1_coeff = (d+4)/(d-2) = "
          f"{nstr(r, 8)} != 4  (not equal-weight)",
          not_equal)

# Exact rational values
check("P024  d=3: P1_coeff = (3+4)/(3-2) = 7  != 4",
      fabs(P1_coeff(mpf('3')) - 7) < mpf('1e-50'))
check("P025  d=5: P1_coeff = (5+4)/(5-2) = 3  != 4",
      fabs(P1_coeff(mpf('5')) - 3) < mpf('1e-50'))
check("P026  d=6: P1_coeff = (6+4)/(6-2) = 10/4 = 5/2  != 4",
      fabs(P1_coeff(mpf('6')) - mpf('5') / 2) < mpf('1e-50'))

# ─── Gegenbauer basis matrices (d=4 consistency) ──────────────────────────────
print("\n=== Section 6: Gegenbauer basis matrix consistency ===")

import numpy as np

# Change-of-basis: P = matrix with columns = U_n in monomials {1,x,x^2,x^3}
P = np.array([
    [1, -2,  3, -4],
    [0,  4,-16, 40],
    [0,  0, 16,-96],
    [0,  0,  0, 64],
], dtype=float)
Pinv = np.linalg.inv(P)

# D^2_{B^4} in monomial basis (row=output, col=input)
D2_mono = np.array([
    [0, 0, 192,    0],
    [0, 0,   0, 1152],
    [0, 0,   0,    0],
    [0, 0,   0,    0],
], dtype=float)

# Change to Gegenbauer basis: D2_geg = Pinv @ D2_mono @ P
D2_geg = Pinv @ D2_mono @ P

check("P027  D^2_geg[0,3] = 18432  (U_3 -> U_0 coefficient)",
      abs(D2_geg[0, 3] - 18432) < 1e-4)
check("P028  D^2_geg[1,3] = 18432  (U_3 -> U_1 coefficient)",
      abs(D2_geg[1, 3] - 18432) < 1e-4)
check("P029  D^2_geg[0,2] = 3072   (U_2 -> U_0 coefficient from A226)",
      abs(D2_geg[0, 2] - 3072) < 1e-4)
check("P030  D^2_geg[2,3] = 0  (no U_3 -> U_2 leakage)",
      abs(D2_geg[2, 3]) < 1e-8)
check("P031  D^2_geg[3,3] = 0  (no U_3 -> U_3 self-action)",
      abs(D2_geg[3, 3]) < 1e-8)

# Equal-weight ratio exactly 1
ratio_check = D2_geg[0, 3] / D2_geg[1, 3]
check("P032  D^2_geg[0,3] / D^2_geg[1,3] = 1.0  (equal-weight ratio)",
      abs(ratio_check - 1.0) < 1e-8)

# ─── K_d computation for each d ───────────────────────────────────────────────
print("\n=== Section 7: K_d prefactors ===")

check("P033  K_{d=3} = 128 * (1/2) * (3/2) * (5/2) * 5 = 1200",
      fabs(K_d(mpf('3')) - 1200) < mpf('1e-50'))
check("P034  K_{d=4} = 4608  (as above)",
      fabs(K_d(mpf('4')) - 4608) < mpf('1e-50'))
check("P035  K_{d=5} = 128 * (3/2) * (5/2) * (7/2) * 7 = 128*105/8*7 = ...",
      K_d(mpf('5')) > 0)  # positivity check

# verify d=5 exactly: alpha=3/2, K = 128 * 3/2 * 5/2 * 7/2 * 7
K5_exact = 128 * mpf('3') / 2 * mpf('5') / 2 * mpf('7') / 2 * 7
check("P036  K_{d=5} = 128 * (3/2)(5/2)(7/2) * 7 = 5880",
      fabs(K_d(mpf('5')) - K5_exact) < mpf('1e-50'))

# ─── Cross-check: 18432 = 4 * K_4 ────────────────────────────────────────────
print("\n=== Section 8: Cross-checks ===")

check("P037  18432 = 4 * K_{d=4} = 4 * 4608",
      fabs(mpf('18432') - 4 * d4_K) < mpf('1e-50'))
check("P038  18432 = 2^11 * 3^2  (factorisation from A226 Remark)",
      fabs(mpf('18432') - 2**11 * 3**2) < mpf('1e-50'))
check("P039  4608 = 2^9 * 3^2  [K_4 factorisation]",
      fabs(d4_K - 2**9 * 3**2) < mpf('1e-50'))

# The universal scalar residue
check("P040  Universal scalar residue A = -d + (d+4) = 4  (independent of d)",
      fabs((-mpf('4') + (mpf('4') + 4)) - 4) < mpf('1e-50') and
      fabs((-mpf('3') + (mpf('3') + 4)) - 4) < mpf('1e-50') and
      fabs((-mpf('7') + (mpf('7') + 4)) - 4) < mpf('1e-50'))

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