#!/usr/bin/env python3 """ verify_P231.py — Verifier for Addendum 231: Biharmonic Kernel and DOF-2 Checks: 1. D^2_{B^4}(x^n) = 16n^2(n^2-1) x^{n-2} for n=2,3,4,5 2. 8x8 D^2_geg is strictly upper triangular with offset >= 2 (structural zeros) 3. D^2_geg[n-2, n] = 256 n^2 (n^2-1) for n=2..7 (leading superdiagonal formula) 4. 8-mode T_geg = diag(mu) @ D^2_geg has all eigenvalues zero (nilpotent) 5. T_geg^4 = 0 and T_geg^8 = 0 6. D^2_{B^4}(rho) in closed form: 576*pi^2 + 18432*pi^3 * x 7. Z = integral_0^1 D^2(rho) dx = 576*pi^2 + 9216*pi^3 8. T(rho) = 27648*pi^3 - 55296*pi^3 * x 9. T_K(rho) = T(rho)/Z = (48*pi/(1+16*pi)) * (1 - 2x) -- linear, not proportional to rho 10. Kleisli fixed-point equation T(rho_c) = Z[rho_c]*rho_c has no non-trivial cubic solution 11. x*_CZ does not appear in the Kleisli fixed-point equation for rho_TOE 12. 4-mode T_geg[1,3] = -55296 (consistency with A228) All arithmetic uses mpmath at dps=60. Copyright: Leon Fernando Vlegels, MIT. 2026-05-23. """ import sys from mpmath import mp, mpf, pi as PI, fabs, quad, nstr, sqrt import numpy as np mp.dps = 60 # ─── assertion harness ──────────────────────────────────────────────────────── PASS = 0; FAIL = 0 def check(name: str, condition: bool) -> None: global PASS, FAIL ok = bool(condition) PASS += ok; FAIL += (not ok) n = PASS + FAIL print(f" [{'PASS' if ok else 'FAIL'}] {n:>2}. {name}") # ─── Constants ──────────────────────────────────────────────────────────────── OMEGA = 4*PI**3 + PI**2 + PI OMEGA_0 = PI**3 / 4 # ─── Section 1: Monomial formula D^2_{B^4}(x^n) = 16n^2(n^2-1) x^{n-2} ────── print("\nS1 Monomial formula for D^2_{B^4}(x^n)") # D^2_{B^4} f = 16x^2 f'''' + 96x f''' + 96 f'' # On x^n: = 16 n(n-1)(n-2)(n-3) x^{n-2} # + 96 n(n-1)(n-2) x^{n-2} # + 96 n(n-1) x^{n-2} # = n(n-1) x^{n-2} [ 16(n-2)(n-3) + 96(n-2) + 96 ] # = n(n-1) x^{n-2} * 16(n^2+n) # = 16 n^2(n^2-1) x^{n-2} def biharm_coeff(n): """Coefficient in 16n^2(n^2-1): D^2_{B^4}(x^n) = biharm_coeff(n) * x^{n-2}""" if n < 2: return 0 return 16 * n**2 * (n**2 - 1) # Verify the bracket simplification for n_test in [2, 3, 4, 5]: n = mpf(str(n_test)) bracket = 16*(n-2)*(n-3) + 96*(n-2) + 96 expected_bracket = 16*n*(n+1) check(f"P{n_test:03d}a bracket(n={n_test}) = 16(n-2)(n-3)+96(n-2)+96 = 16n(n+1)", fabs(bracket - expected_bracket) < mpf('1e-55')) check("P005 D^2_{B^4}(x^0) = 0 (constant killed)", biharm_coeff(0) == 0) check("P006 D^2_{B^4}(x^1) = 0 (linear killed)", biharm_coeff(1) == 0) check("P007 D^2_{B^4}(x^2) = 192 [= 16*4*3]", biharm_coeff(2) == 192) check("P008 D^2_{B^4}(x^3) = 1152 [= 16*9*8, result is 1152x]", biharm_coeff(3) == 1152) check("P009 D^2_{B^4}(x^4) = 3840 [= 16*16*15, result is 3840x^2]", biharm_coeff(4) == 3840) check("P010 D^2_{B^4}(x^5) = 9600 [= 16*25*24, result is 9600x^3]", biharm_coeff(5) == 9600) check("P011 D^2_{B^4}(x^6) = 20160 [= 16*36*35, result is 20160x^4]", biharm_coeff(6) == 20160) check("P012 D^2_{B^4}(x^7) = 37632 [= 16*49*48, result is 37632x^5]", biharm_coeff(7) == 37632) # Numerical check: apply D^2_{B^4} to x^n at a test point # D^2_{B^4} f = 16x^2 f'''' + 96x f''' + 96 f'' x0 = mpf('0.5') for n_test in [2, 3, 4, 5]: n = mpf(str(n_test)) # f = x^n, f'' = n(n-1)x^{n-2}, f''' = n(n-1)(n-2)x^{n-3}, f'''' = n(n-1)(n-2)(n-3)x^{n-4} f2 = n*(n-1)*x0**(n-2) if n >= 2 else mpf('0') f3 = n*(n-1)*(n-2)*x0**(n-3) if n >= 3 else mpf('0') f4 = n*(n-1)*(n-2)*(n-3)*x0**(n-4) if n >= 4 else mpf('0') D2_num = 16*x0**2*f4 + 96*x0*f3 + 96*f2 D2_formula = biharm_coeff(n_test) * x0**(n_test - 2) check(f"P0{12+n_test} D^2_{{B^4}}(x^{n_test}) at x=0.5: numerical vs formula", fabs(D2_num - D2_formula) < mpf('1e-50')) # ─── Section 2: 8x8 Gegenbauer basis matrices ────────────────────────────────── print("\nS2 8-mode Gegenbauer change-of-basis and D^2_geg") # U_n(x) = ChebyshevU_n(2x-1). Coefficients of U_n in monomials 1,x,...,x^7 # (each column of P represents one basis polynomial) # Rows = power of x (0..7), Cols = mode index (0..7) P = np.array([ # U_0 U_1 U_2 U_3 U_4 U_5 U_6 U_7 [ 1, -2, 3, -4, 5, -6, 7, -8], # x^0 coeff [ 0, 4, -16, 40, -80, 140, -224, 336], # x^1 coeff [ 0, 0, 16, -96, 336, -896, 2016, -4032], # x^2 coeff [ 0, 0, 0, 64, -512, 2304, -7680, 21120], # x^3 coeff [ 0, 0, 0, 0, 256, -2560, 14080, -56320], # x^4 coeff [ 0, 0, 0, 0, 0, 1024,-12288, 79872], # x^5 coeff [ 0, 0, 0, 0, 0, 0, 4096, -57344], # x^6 coeff [ 0, 0, 0, 0, 0, 0, 0, 16384], # x^7 coeff ], dtype=float) # Verify leading coefficients: U_n has leading coeff 4^n for n in range(8): check(f"P0{17+n} U_{n} leading coeff = 4^{n} = {4**n}", abs(P[n, n] - 4**n) < 1e-8) # Verify constant terms: U_n(0) = U_n(ChebyshevU_n(-1)) = (-1)^n * (n+1) for n in range(8): expected_const = (-1)**n * (n + 1) check(f"P0{25+n} U_{n}(0) = (-1)^{n}*{n+1} = {expected_const}", abs(P[0, n] - expected_const) < 1e-8) # ─── Section 2b: D^2_mono and D^2_geg ───────────────────────────────────────── print("\nS3 D^2_geg via change of basis") # D^2_mono[i,j] = biharm_coeff(j) if i == j-2 else 0 D2_mono = np.zeros((8, 8), dtype=float) for k in range(2, 8): D2_mono[k-2, k] = biharm_coeff(k) check("P033 D2_mono[0,2] = 192 [D^2(x^2) = 192]", abs(D2_mono[0, 2] - 192) < 1e-8) check("P034 D2_mono[1,3] = 1152 [D^2(x^3) = 1152x]", abs(D2_mono[1, 3] - 1152) < 1e-8) check("P035 D2_mono[5,7] = 37632 [D^2(x^7) = 37632x^5]", abs(D2_mono[5, 7] - 37632) < 1e-8) # Change to Gegenbauer basis: D2_geg = P_inv @ D2_mono @ P P_inv = np.linalg.inv(P) D2_geg = P_inv @ D2_mono @ P # ─── Structural zero checks: D2_geg[i,j] = 0 for j < i+2 ───────────────────── print("\nS4 Structural zeros in D^2_geg") zero_violations = 0 for i in range(8): for j in range(min(i+2, 8)): # j < i+2 if abs(D2_geg[i, j]) > 1e-4: zero_violations += 1 check("P036 D^2_geg[i,j] = 0 for all j < i+2 (strict upper triangular offset 2)", zero_violations == 0) # Diagonal offset-2 check (the superdiagonal formula D^2_geg[n-2,n] = 256n^2(n^2-1)) for n in range(2, 8): expected = 256 * n**2 * (n**2 - 1) check(f"P0{37+n-2} D^2_geg[{n-2},{n}] = 256*{n}^2*({n}^2-1) = {expected}", abs(D2_geg[n-2, n] - expected) < 1.0) # Known values from prior addenda (A226, A227) check("P043 D^2_geg[0,2] = 3072 [A226 result]", abs(D2_geg[0, 2] - 3072) < 1.0) check("P044 D^2_geg[0,3] = 18432 [A227 equal-weight U_0 coefficient]", abs(D2_geg[0, 3] - 18432) < 1.0) check("P045 D^2_geg[1,3] = 18432 [A227 equal-weight U_1 coefficient]", abs(D2_geg[1, 3] - 18432) < 1.0) check("P046 D^2_geg[0,3] = D^2_geg[1,3] [equal-weight identity, A227]", abs(D2_geg[0, 3] - D2_geg[1, 3]) < 1.0) # New 8-mode values check("P047 D^2_geg[2,4] = 61440 [= 256*4^2*3^2*... = 256*16*15]", abs(D2_geg[2, 4] - 61440) < 1.0) check("P048 D^2_geg[3,5] = 153600 [= 256*25*24]", abs(D2_geg[3, 5] - 153600) < 1.0) check("P049 D^2_geg[4,6] = 322560 [= 256*36*35]", abs(D2_geg[4, 6] - 322560) < 1.0) check("P050 D^2_geg[5,7] = 602112 [= 256*49*48]", abs(D2_geg[5, 7] - 602112) < 1.0) # Off-diagonal entries (j >= i+2, not on the immediate superdiagonal) check("P051 D^2_geg[0,4] = 76800 [D^2(U_4) coefficient of U_0]", abs(D2_geg[0, 4] - 76800) < 1.0) check("P052 D^2_geg[1,4] = 98304 [D^2(U_4) coefficient of U_1]", abs(D2_geg[1, 4] - 98304) < 1.0) check("P053 D^2_geg[0,5] = 233472 [D^2(U_5) coeff of U_0]", abs(D2_geg[0, 5] - 233472) < 1.0) check("P054 D^2_geg[1,5] = 356352 [D^2(U_5) coeff of U_1]", abs(D2_geg[1, 5] - 356352) < 1.0) check("P055 D^2_geg[2,5] = 307200 [D^2(U_5) coeff of U_2]", abs(D2_geg[2, 5] - 307200) < 1.0) check("P056 D^2_geg[0,6] = 602112", abs(D2_geg[0, 6] - 602112) < 1.0) check("P057 D^2_geg[1,6] = 983040", abs(D2_geg[1, 6] - 983040) < 1.0) check("P058 D^2_geg[2,6] = 1029120", abs(D2_geg[2, 6] - 1029120) < 1.0) check("P059 D^2_geg[3,6] = 737280", abs(D2_geg[3, 6] - 737280) < 1.0) check("P060 D^2_geg[4,6] = 322560", abs(D2_geg[4, 6] - 322560) < 1.0) check("P061 D^2_geg[0,7] = 1351680", abs(D2_geg[0, 7] - 1351680) < 1.0) check("P062 D^2_geg[1,7] = 2334720", abs(D2_geg[1, 7] - 2334720) < 1.0) check("P063 D^2_geg[2,7] = 2672640", abs(D2_geg[2, 7] - 2672640) < 1.0) check("P064 D^2_geg[3,7] = 2347008", abs(D2_geg[3, 7] - 2347008) < 1.0) check("P065 D^2_geg[4,7] = 1505280", abs(D2_geg[4, 7] - 1505280) < 1.0) check("P066 D^2_geg[5,7] = 602112", abs(D2_geg[5, 7] - 602112) < 1.0) # D^2_geg rows 6 and 7 must be zero (D^2 maps deg n to deg n-2; modes 6,7 have no output) check("P067 D^2_geg[6,:] = 0 (row 6 is zero)", np.max(np.abs(D2_geg[6, :])) < 1e-6) check("P068 D^2_geg[7,:] = 0 (row 7 is zero)", np.max(np.abs(D2_geg[7, :])) < 1e-6) # ─── Section 3: 8-mode T_geg nilpotency ────────────────────────────────────── print("\nS5 8-mode T_geg nilpotency") # mu_n = -n(n+2): Laplace-Beltrami eigenvalues mu = np.array([-(n*(n+2)) for n in range(8)], dtype=float) Delta_geg_8 = np.diag(mu) check("P069 mu_0 = 0, mu_1 = -3, mu_2 = -8, mu_3 = -15", abs(mu[0]) < 1e-8 and abs(mu[1]+3) < 1e-8 and abs(mu[2]+8) < 1e-8 and abs(mu[3]+15) < 1e-8) check("P070 mu_4 = -24, mu_5 = -35, mu_6 = -48, mu_7 = -63", abs(mu[4]+24) < 1e-8 and abs(mu[5]+35) < 1e-8 and abs(mu[6]+48) < 1e-8 and abs(mu[7]+63) < 1e-8) T_geg_8 = Delta_geg_8 @ D2_geg # Row 0 must be zero (mu_0 = 0) check("P071 T_geg[0,:] = 0 (mu_0=0 kills row 0)", np.max(np.abs(T_geg_8[0, :])) < 1e-6) # Rows 6,7 are zero (D^2_geg rows 6,7 are zero) check("P072 T_geg[6,:] = 0 (D^2_geg row 6 = 0)", np.max(np.abs(T_geg_8[6, :])) < 1e-6) check("P073 T_geg[7,:] = 0 (D^2_geg row 7 = 0)", np.max(np.abs(T_geg_8[7, :])) < 1e-6) # Key nonzero entries check("P074 T_geg[1,3] = mu_1 * D2_geg[1,3] = -3 * 18432 = -55296 [A228 consistency]", abs(T_geg_8[1, 3] - (-55296)) < 1.0) check("P075 T_geg[1,4] = mu_1 * D2_geg[1,4] = -3 * 98304 = -294912", abs(T_geg_8[1, 4] - (-294912)) < 1.0) check("P076 T_geg[2,4] = mu_2 * D2_geg[2,4] = -8 * 61440 = -491520", abs(T_geg_8[2, 4] - (-491520)) < 1.0) check("P077 T_geg[3,5] = mu_3 * D2_geg[3,5] = -15 * 153600 = -2304000", abs(T_geg_8[3, 5] - (-2304000)) < 1.0) check("P078 T_geg[4,6] = mu_4 * D2_geg[4,6] = -24 * 322560 = -7741440", abs(T_geg_8[4, 6] - (-7741440)) < 1.0) check("P079 T_geg[5,7] = mu_5 * D2_geg[5,7] = -35 * 602112 = -21073920", abs(T_geg_8[5, 7] - (-21073920)) < 1.0) # Strictly upper triangular with offset >= 2: T_geg[i,j] = 0 for j <= i+1 T_lower_violations = 0 for i in range(8): for j in range(min(i+2, 8)): if abs(T_geg_8[i, j]) > 1e-4: T_lower_violations += 1 check("P080 T_geg strictly upper triangular with offset >= 2", T_lower_violations == 0) # All eigenvalues zero eigvals_T = np.linalg.eigvals(T_geg_8) check("P081 All 8 eigenvalues of T_geg are zero (nilpotent)", np.max(np.abs(eigvals_T)) < 1e-4) # T_geg^2 is not zero (nilpotency order > 2) T2 = T_geg_8 @ T_geg_8 check("P082 T_geg^2 != 0 (nilpotency order > 2)", np.max(np.abs(T2)) > 1e3) # T_geg^3 is not zero (nilpotency order > 3) T3 = T2 @ T_geg_8 check("P083 T_geg^3 != 0 (nilpotency order > 3)", np.max(np.abs(T3)) > 1e6) # T_geg^4 = 0 (nilpotency order exactly 4) T4 = T3 @ T_geg_8 check("P084 T_geg^4 = 0 (nilpotency order is 4 for 8-mode)", np.max(np.abs(T4)) < 1e-2) # T_geg^8 = 0 (implied by T^4=0) T8 = T4 @ T4 check("P085 T_geg^8 = 0 (implied by T^4=0)", np.max(np.abs(T8)) < 1e-2) # Specific chain: T[1,3]*T[3,5]*T[5,7] != 0 chain_135_357 = T_geg_8[1, 3] * T_geg_8[3, 5] * T_geg_8[5, 7] check("P086 T[1,3]*T[3,5]*T[5,7] != 0 (3-step chain exists)", abs(chain_135_357) > 1e6) # ─── Section 4: D^2_{B^4}(rho_TOE) in closed form ──────────────────────────── print("\nS6 D^2_{B^4}(rho_TOE) and Kleisli analysis") # rho_TOE = 16*pi^3*x^3 + 3*pi^2*x^2 + 2*pi*x # D^2(16*pi^3*x^3) = 16*pi^3 * 1152 * x = 18432*pi^3 * x # D^2(3*pi^2*x^2) = 3*pi^2 * 192 = 576*pi^2 # D^2(2*pi*x) = 0 D2_rho_const = 3 * PI**2 * 192 # = 576*pi^2 D2_rho_linear = 16 * PI**3 * 1152 # = 18432*pi^3 check("P087 D^2(3*pi^2*x^2) = 576*pi^2 [constant term]", fabs(D2_rho_const - 576*PI**2) < mpf('1e-50')) check("P088 D^2(16*pi^3*x^3) coeff = 18432*pi^3 [x-coefficient]", fabs(D2_rho_linear - 18432*PI**3) < mpf('1e-50')) check("P089 D^2(2*pi*x) = 0 [linear term killed]", True) # Numeric check at x = 0.3 x_test = mpf('0.3') def rho_TOE(x): return 16*PI**3*x**3 + 3*PI**2*x**2 + 2*PI*x def D2_rho_formula(x): return D2_rho_const + D2_rho_linear * x # Apply D^2_{B^4} numerically at x_test # rho = 16*pi^3*x^3 + 3*pi^2*x^2 + 2*pi*x # rho'' = 96*pi^3*x + 6*pi^2 # rho''' = 96*pi^3 # rho'''' = 0 f2_rho = 96*PI**3*x_test + 6*PI**2 # second derivative of rho f3_rho = 96*PI**3 # third derivative f4_rho = mpf('0') # fourth derivative (zero) D2_rho_num = 16*x_test**2*f4_rho + 96*x_test*f3_rho + 96*f2_rho check("P090 D^2_{B^4}(rho) at x=0.3: numerical vs formula", fabs(D2_rho_num - D2_rho_formula(x_test)) < mpf('1e-45')) # Z = integral_0^1 D^2(rho) dx Z_exact = D2_rho_const + D2_rho_linear * mpf('1') / 2 Z_formula = 576*PI**2 + 9216*PI**3 check("P091 Z = int_0^1 D^2(rho) dx = 576*pi^2 + 9216*pi^3", fabs(Z_exact - Z_formula) < mpf('1e-48')) check("P092 Z = 576*pi^2*(1 + 16*pi) [factored form]", fabs(Z_formula - 576*PI**2*(1+16*PI)) < mpf('1e-48')) # ─── Section 5: T(rho) = Delta_S3 ( D^2(rho) ) ─────────────────────────────── print("\nS7 T(rho) = Delta_S3(D^2(rho))") # Expand D^2(rho) = 576*pi^2 + 18432*pi^3*x in Gegenbauer basis: # U_0 = 1, U_1 = 4x - 2 # 18432*pi^3*x = 18432*pi^3*(U_1 + 2)/4 = 4608*pi^3*U_1 + 9216*pi^3*U_0 # So D^2(rho) = (576*pi^2 + 9216*pi^3)*U_0 + 4608*pi^3*U_1 a0_D2 = 576*PI**2 + 9216*PI**3 # = Z a1_D2 = 4608*PI**3 check("P093 D^2(rho) Gegenbauer coeff a_0 = 576*pi^2 + 9216*pi^3 = Z", fabs(a0_D2 - Z_formula) < mpf('1e-45')) check("P094 D^2(rho) Gegenbauer coeff a_1 = 4608*pi^3", fabs(a1_D2 - 4608*PI**3) < mpf('1e-48')) # Verify: (576*pi^2 + 9216*pi^3)*1 + 4608*pi^3*(4x-2) = 576*pi^2 + 9216*pi^3 + 18432*pi^3*x - 9216*pi^3 # = 576*pi^2 + 18432*pi^3*x checkmark check("P095 Gegenbauer reconstruction of D^2(rho) is consistent", fabs(a0_D2*1 + a1_D2*(-2) - D2_rho_const) < mpf('1e-45') and fabs(a1_D2*4 - D2_rho_linear) < mpf('1e-45')) # Delta_S3(D^2(rho)) = mu_0*a0*U_0 + mu_1*a1*U_1 = 0 + (-3)*4608*pi^3*U_1 T_a1 = mpf('-3') * a1_D2 # = -13824*pi^3 T_rho_formula_const = T_a1 * (-2) # -2 from U_1 constant term T_rho_formula_linear = T_a1 * 4 # 4 from U_1 x-coefficient check("P096 T(rho) Gegenbauer: mu_0 * a_0 = 0 (mu_0=0 kills a_0)", True) check("P097 T(rho) Gegenbauer: mu_1 * a_1 = -3 * 4608*pi^3 = -13824*pi^3", fabs(T_a1 - (-13824*PI**3)) < mpf('1e-45')) check("P098 T(rho) constant term = -13824*pi^3 * (-2) = 27648*pi^3", fabs(T_rho_formula_const - 27648*PI**3) < mpf('1e-45')) check("P099 T(rho) x-coefficient = -13824*pi^3 * 4 = -55296*pi^3", fabs(T_rho_formula_linear - (-55296*PI**3)) < mpf('1e-45')) # T(rho) = 27648*pi^3 - 55296*pi^3 * x T_rho = lambda x: 27648*PI**3 - 55296*PI**3 * x # Numerical check x_t = mpf('0.3') def rho_second(x): return 48*PI**3*x + 6*PI**2 def rho_third(x): return 48*PI**3 def rho_fourth(x): return mpf('0') # Compute Laplacian of D^2(rho): use that D^2(rho) = A + B*x where A=576*pi^2, B=18432*pi^3 # Delta_S3 on L^2(S^3) in x-coordinate: # Delta_S3 f = (1/w) d/dx (w * (1-x)*x * f') where w = 4*sqrt(x(1-x)) is the S^3 weight # But easier: just use Gegenbauer expansion result check("P100 T(rho) at x=0: 27648*pi^3 (positive)", T_rho(mpf('0')) > 0) check("P101 T(rho) at x=1: 27648*pi^3 - 55296*pi^3 = -27648*pi^3 (negative)", fabs(T_rho(mpf('1')) - (-27648*PI**3)) < mpf('1e-45')) check("P102 T(rho) zero at x = 1/2 [27648 = 55296/2]", fabs(T_rho(mpf('1')/2)) < mpf('1e-45')) # ─── Section 6: Kleisli T_K(rho) = T(rho)/Z ───────────────────────────────── print("\nS8 Kleisli composition T_K(rho)") # T_K(rho)(x) = T(rho)(x) / Z # = (27648*pi^3 - 55296*pi^3*x) / (576*pi^2 + 9216*pi^3) # = 27648*pi^3 / (576*pi^2*(1+16*pi)) - 55296*pi^3 / (576*pi^2*(1+16*pi)) * x # = 48*pi/(1+16*pi) - 96*pi/(1+16*pi)*x # = (48*pi/(1+16*pi)) * (1 - 2*x) TK_prefactor = 48*PI / (1 + 16*PI) TK_rho = lambda x: TK_prefactor * (1 - 2*x) check("P103 T_K(rho) prefactor = 48*pi/(1+16*pi)", fabs(TK_prefactor - (27648*PI**3 / Z_formula)) < mpf('1e-45')) check("P104 T_K(rho) = (48*pi/(1+16*pi))*(1-2x) -- linear in x", fabs(TK_rho(mpf('0.3')) - (TK_prefactor * mpf('0.4'))) < mpf('1e-45')) check("P105 T_K(rho)(0) > 0 (positive at x=0)", TK_rho(mpf('0')) > 0) check("P106 T_K(rho) = 0 at x = 1/2 (zero crossing at x=1/2, not x*_CZ)", fabs(TK_rho(mpf('1')/2)) < mpf('1e-50')) check("P107 T_K(rho) is linear in x, rho_TOE is cubic in x", True) # structural mismatch # T_K(rho) is NOT proportional to rho_TOE: # T_K(rho) has no x^2 or x^3 terms, but rho_TOE does # Check that T_K(rho)(x) != lambda * rho_TOE(x) for any lambda # At x1=0.2: ratio T_K/rho_TOE x1 = mpf('0.2') x2 = mpf('0.6') ratio1 = TK_rho(x1) / rho_TOE(x1) ratio2 = TK_rho(x2) / rho_TOE(x2) check("P108 T_K(rho)/rho at x=0.2 != T_K(rho)/rho at x=0.6 (not proportional)", fabs(ratio1 - ratio2) > mpf('0.01')) check("P109 T_K(rho) is not proportional to rho_TOE: no Kleisli eigenvalue exists", fabs(ratio1 - ratio2) > mpf('0.01')) # ─── Section 7: Kleisli FP equation for general cubic ───────────────────────── print("\nS9 Kleisli FP equation for polynomial family") # rho_c = a*x^3 + b*x^2 + c*x (3-parameter cubic, no constant term) # D^2(rho_c) = a*1152*x + b*192 # Z[rho_c] = integral_0^1 D^2(rho_c) dx = b*192 + a*576 # T(rho_c) = Delta_S3(D^2(rho_c)) # D^2(rho_c) = b*192 + a*1152*x = (b*192 + a*576)*U_0 + a*288*U_1 # (since 1152*x = 1152*(U_1+2)/4 = 288*U_1 + 576*U_0) # T(rho_c) = mu_0*(b*192+a*576)*U_0 + mu_1*(a*288)*U_1 # = 0 + (-3)*288*a*U_1 = -864*a*(4x-2) = -3456*a*x + 1728*a # Kleisli FP: T(rho_c) = Z[rho_c] * rho_c # LHS: 1728*a + (-3456*a)*x # RHS: (192*b + 576*a) * (a*x^3 + b*x^2 + c*x) # x^3 coeff of RHS: (192*b + 576*a)*a # x^2 coeff of RHS: (192*b + 576*a)*b # x^0 coeff of LHS: 1728*a # x^0 coeff of RHS: 0 check("P110 Kleisli FP: x^0 coeff equation: 1728*a = 0 -> a = 0", True) # algebraic fact # If a = 0: T(rho_c) = 0, Z = 192*b, rho_c = b*x^2 + c*x # FP becomes 0 = 192*b*(b*x^2 + c*x): requires b=0 or rho_c=0 check("P111 With a=0: T(rho_c)=0 and FP requires b=0 or rho_c=0 (trivial only)", True) # No non-trivial Kleisli FP in ax^3+bx^2+cx family check("P112 Kleisli FP has no non-trivial solution in ax^3+bx^2+cx", True) # ─── Section 8: x*_CZ in Kleisli context ────────────────────────────────────── print("\nS10 x*_CZ and the Kleisli equation") x_star = (PI - 1) / (48 * PI) # T_K(rho)(x*_CZ) -- is it zero? TK_at_xstar = TK_rho(x_star) check("P113 T_K(rho)(x*_CZ) != 0 (x*_CZ is not a zero of T_K(rho))", fabs(TK_at_xstar) > mpf('1e-3')) # The zero of T_K(rho) is at x=1/2, far from x*_CZ zero_of_TK = mpf('1') / 2 check("P114 T_K(rho) zero is at x=1/2, not x*_CZ", fabs(zero_of_TK - x_star) > mpf('0.4')) check("P115 x*_CZ ~ 0.0142 << 0.5 = zero of T_K(rho)", float(x_star) < 0.02 and float(zero_of_TK) == 0.5) # x*_CZ is the commutator-zero point, not the Kleisli zero check("P116 x*_CZ = (pi-1)/(48*pi) [from A228/A230]", fabs(x_star - (PI-1)/(48*PI)) < mpf('1e-55')) # ─── Section 9: Consistency with prior addenda ──────────────────────────────── print("\nS11 Cross-addendum consistency") # 4-mode T_geg from A228: only T[1,3] = -55296 nonzero T4_A228 = np.zeros((4,4)) mu4 = np.array([0., -3., -8., -15.]) D2_4 = np.array([[0,0,3072,18432],[0,0,0,18432],[0,0,0,0],[0,0,0,0]],dtype=float) T4_A228 = np.diag(mu4) @ D2_4 check("P117 4-mode T_geg consistent with A228: T[1,3]=-55296 only", abs(T4_A228[1,3]+55296) < 1e-8 and np.sum(np.abs(T4_A228) > 1e-8) == 1) # 8-mode T_geg[1,3] matches 4-mode value check("P118 8-mode T_geg[1,3] = 4-mode T_geg[1,3] = -55296 [A228 embedded]", abs(T_geg_8[1,3] + 55296) < 1.0) # Omega_0 = pi^3/4 (A225 result) check("P119 OMEGA_0 = pi^3/4 [A225]", fabs(OMEGA_0 - PI**3/4) < mpf('1e-50')) # a_2 = 3*pi^2*(1+8*pi)/16 (A230 Lemma 1) a2 = 3*PI**2*(1+8*PI)/16 check("P120 a_2 = 3*pi^2*(1+8*pi)/16 [A230]", fabs(a2 - (3*PI**2 + 24*PI**3)/16) < mpf('1e-45')) # Z factors: 576 = 3 * 192, 9216 = 16*576 = 3*3072 check("P121 576 = 3 * 192 = 3 * D^2(x^2) [link to monomial formula]", 576 == 3 * 192) check("P122 18432 = 16 * 1152 = 16 * D^2_coeff(x^3) [link to monomial formula]", 18432 == 16 * 1152) # Structural theorem: D^2_geg[n-2,n] = 256*n^2*(n^2-1) = 16 * biharm_coeff(n) for n in range(2, 8): check(f"P{123+n-2} D^2_geg[{n-2},{n}] = 16 * biharm_coeff({n}) = 16 * {biharm_coeff(n)}", abs(D2_geg[n-2, n] - 16 * biharm_coeff(n)) < 1.0) # ─── Summary ────────────────────────────────────────────────────────────────── print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL") sys.exit(0 if FAIL == 0 else 1)