""" verify_P226.py — Verification suite for Addendum 226 "Wheel Commutator [D²_B⁴, Δ_S³] in the Gegenbauer Basis" L. F. Vlegels, 22 May 2026 © Léon Fernando Vlegels. MIT License. Assertions: P001–P008 D²_B⁴ matrix in {1, x, x², x³} basis P009–P016 Δ_S³ matrix in {1, x, x², x³} basis: values and eigenvalue check P017–P024 Commutator [D²_B⁴, Δ_S³] in monomial basis P025–P030 D²_B⁴ in Gegenbauer basis: equal-weight mixing identity P031–P037 Spectral-gap formula: commutator = D²·(μⱼ − μᵢ) P038–P042 Layer descent: [D²,Δ]ρ ∈ span{1,x} P043–P046 Robertson bound: ratio far from unity (HUP not saturated) Expected outcome: 46 assertions, all passing. """ import sys import os import math import numpy as np from fractions import Fraction # ── mpmath ──────────────────────────────────────────────────────────────────── from mpmath import mp, mpf, pi as mpi, fabs, nstr mp.dps = 60 # ── Path setup ──────────────────────────────────────────────────────────────── REPO_ROOT = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", "..", "..")) if REPO_ROOT not in sys.path: sys.path.insert(0, REPO_ROOT) # ── Assertion infrastructure ────────────────────────────────────────────────── PASS = 0 FAIL = 0 _N = 0 def check(name: str, condition: bool) -> None: global PASS, FAIL, _N _N += 1 ok = bool(condition) PASS += ok FAIL += (not ok) print(f" [{'PASS' if ok else 'FAIL'}] {_N:>2}. {name}") # ── Constants ───────────────────────────────────────────────────────────────── PI = float(mpi) OMEGA = 4*PI**3 + PI**2 + PI OMEGA_0 = PI**3 / 4 # kernel constant # Operator matrices in {1, x, x², x³} (exact, as numpy int arrays) D2B4_mono = np.array([ [0, 0, 192, 0], [0, 0, 0, 1152], [0, 0, 0, 0], [0, 0, 0, 0], ], dtype=float) DS3_mono = np.array([ [0, 1.5, 0, 0], [0, -3, 5, 0], [0, 0, -8, 10.5], [0, 0, 0, -15], ], dtype=float) # Change-of-basis: P[:,n] = U_n in monomial coords 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) # Operators in Gegenbauer basis D2B4_geg = Pinv @ D2B4_mono @ P DS3_geg = Pinv @ DS3_mono @ P COM_mono = D2B4_mono @ DS3_mono - DS3_mono @ D2B4_mono COM_geg = D2B4_geg @ DS3_geg - DS3_geg @ D2B4_geg # ρ coefficient vector in monomial basis: [0, 2π, 3π², 16π³] rho_mono = np.array([0, 2*PI, 3*PI**2, 16*PI**3]) COM_rho = COM_mono @ rho_mono # [D²,Δ]ρ in monomial basis # ───────────────────────────────────────────────────────────────────────────── # P001–P008 D²_B⁴ matrix in monomial basis # ───────────────────────────────────────────────────────────────────────────── print("\n=== P001–P008 D²_B⁴ matrix in {1, x, x², x³} ===") # Formula: D²_B⁴ f(x) = 16x² f'''' + 96x f''' + 96 f'' check("P001 D²_B⁴(1) = 0 (constant is annihilated)", abs(D2B4_mono[0,0]) < 1e-12 and abs(D2B4_mono[1,0]) < 1e-12) check("P002 D²_B⁴(x) = 0 (linear is annihilated)", abs(D2B4_mono[0,1]) < 1e-12 and abs(D2B4_mono[1,1]) < 1e-12) check("P003 D²_B⁴(x²) = 192 (constant output)", abs(D2B4_mono[0,2] - 192) < 1e-10) check("P004 D²_B⁴(x²) has no x component (output is pure constant)", abs(D2B4_mono[1,2]) < 1e-12) check("P005 D²_B⁴(x³) = 1152·x (linear output)", abs(D2B4_mono[1,3] - 1152) < 1e-10) check("P006 D²_B⁴(x³) has no constant component", abs(D2B4_mono[0,3]) < 1e-12) check("P007 D²_B⁴ is strictly lower-triangular in monomial basis (rows 2,3 are zero)", np.all(np.abs(D2B4_mono[2:,:]) < 1e-12)) check("P008 D²_B⁴ annihilates polynomials of degree ≤ 1", np.all(np.abs(D2B4_mono[:,0:2]) < 1e-12)) # ───────────────────────────────────────────────────────────────────────────── # P009–P016 Δ_S³ matrix in monomial basis # ───────────────────────────────────────────────────────────────────────────── print("\n=== P009–P016 Δ_S³ matrix in {1, x, x², x³} ===") check("P009 Δ_S³(1) = 0 (constants are in the kernel)", np.all(np.abs(DS3_mono[:,0]) < 1e-12)) check("P010 Δ_S³(x) = 3/2 - 3x (coeff of 1: 3/2)", abs(DS3_mono[0,1] - 1.5) < 1e-12) check("P011 Δ_S³(x) = 3/2 - 3x (coeff of x: -3)", abs(DS3_mono[1,1] - (-3)) < 1e-12) check("P012 Δ_S³(x²) = 5x - 8x²", abs(DS3_mono[1,2] - 5) < 1e-12 and abs(DS3_mono[2,2] - (-8)) < 1e-12) check("P013 Δ_S³(x³) = (21/2)x² - 15x³", abs(DS3_mono[2,3] - 10.5) < 1e-12 and abs(DS3_mono[3,3] - (-15)) < 1e-12) # Eigenvalue verification: Δ_S³(U_n) = μ_n U_n # U_1 = 4x-2 → Δ_S³(4x-2) = 4(3/2-3x) = 6-12x = -3(4x-2) ✓ U1_coeffs = np.array([-2, 4, 0, 0], dtype=float) DS3_U1 = DS3_mono @ U1_coeffs check("P014 Δ_S³(U_1) = -3·U_1 (eigenvalue μ₁ = -3)", abs(DS3_U1[0] - (-3)*(-2)) < 1e-10 and abs(DS3_U1[1] - (-3)*(4)) < 1e-10) # U_2 = 16x²-16x+3 → eigenvalue -8 U2_coeffs = np.array([3, -16, 16, 0], dtype=float) DS3_U2 = DS3_mono @ U2_coeffs expected_U2 = -8 * U2_coeffs check("P015 Δ_S³(U_2) = -8·U_2 (eigenvalue μ₂ = -8)", np.all(np.abs(DS3_U2 - expected_U2) < 1e-10)) # U_3 = 64x³-96x²+40x-4 → eigenvalue -15 U3_coeffs = np.array([-4, 40, -96, 64], dtype=float) DS3_U3 = DS3_mono @ U3_coeffs expected_U3 = -15 * U3_coeffs check("P016 Δ_S³(U_3) = -15·U_3 (eigenvalue μ₃ = -15)", np.all(np.abs(DS3_U3 - expected_U3) < 1e-8)) # ───────────────────────────────────────────────────────────────────────────── # P017–P024 Commutator in monomial basis # ───────────────────────────────────────────────────────────────────────────── print("\n=== P017–P024 Commutator [D²_B⁴, Δ_S³] in monomial basis ===") check("P017 [D²,Δ](1) = 0 (commutator annihilates constants)", np.all(np.abs(COM_mono[:,0]) < 1e-10)) check("P018 [D²,Δ](x) = 0 (commutator annihilates linear functions)", np.all(np.abs(COM_mono[:,1]) < 1e-10)) check("P019 [D²,Δ](x²): constant entry = -1536", abs(COM_mono[0,2] - (-1536)) < 1e-8) check("P020 [D²,Δ](x²): no x component", abs(COM_mono[1,2]) < 1e-10) check("P021 [D²,Δ](x³): constant entry = 288", abs(COM_mono[0,3] - 288) < 1e-8) check("P022 [D²,Δ](x³): x entry = -13824", abs(COM_mono[1,3] - (-13824)) < 1e-8) check("P023 [D²,Δ] rows 2,3 are identically zero (layer descent)", np.all(np.abs(COM_mono[2:,:]) < 1e-10)) check("P024 All commutator entries divisible by 96 (common factor)", all(abs(v % 96) < 1e-8 for v in [-1536, 288, -13824])) # ───────────────────────────────────────────────────────────────────────────── # P025–P030 D²_B⁴ in Gegenbauer basis: equal-weight mixing identity # ───────────────────────────────────────────────────────────────────────────── print("\n=== P025–P030 D²_B⁴ in Gegenbauer basis ===") # Δ_S³ in Gegenbauer basis should be diagonal(0,-3,-8,-15) check("P025 Δ_S³ in Gegenbauer basis is diagonal (eigenbasis check)", np.all(np.abs(DS3_geg - np.diag([0,-3,-8,-15])) < 1e-8)) check("P026 D²_B⁴(U_0) = 0 (constant mode annihilated)", np.all(np.abs(D2B4_geg[:,0]) < 1e-8)) check("P027 D²_B⁴(U_1) = 0 (dipole mode annihilated)", np.all(np.abs(D2B4_geg[:,1]) < 1e-8)) check("P028 D²_B⁴(U_2) = 3072·U_0 (quadrupole → scalar)", abs(D2B4_geg[0,2] - 3072) < 1e-6 and np.all(np.abs(D2B4_geg[1:,2]) < 1e-6)) check("P029 D²_B⁴(U_3) = 18432·U_0 + 18432·U_1 (equal-weight mixing)", abs(D2B4_geg[0,3] - 18432) < 1e-4 and abs(D2B4_geg[1,3] - 18432) < 1e-4) check("P030 Equal-weight ratio: D²_B⁴(U_3)|_{U_0} / D²_B⁴(U_3)|_{U_1} = 1", abs(D2B4_geg[0,3] / D2B4_geg[1,3] - 1.0) < 1e-8) # ───────────────────────────────────────────────────────────────────────────── # P031–P037 Spectral-gap formula # ───────────────────────────────────────────────────────────────────────────── print("\n=== P031–P037 Spectral-gap formula [D²,Δ]_{ij} = D²_{ij}(μⱼ - μᵢ) ===") mu = lambda n: -n*(n+2) # Verify each non-zero commutator entry # [0,2]: D²_geg[0,2] * (μ_2 - μ_0) = 3072 * (-8 - 0) = -24576 pred_02 = D2B4_geg[0,2] * (mu(2) - mu(0)) check("P031 [D²,Δ]_geg[0,2] = D²[0,2]·(μ₂-μ₀) = 3072·(-8) = -24576", abs(COM_geg[0,2] - pred_02) < 1e-6 and abs(pred_02 - (-24576)) < 1e-6) # [0,3]: D²_geg[0,3] * (μ_3 - μ_0) = 18432 * (-15) = -276480 pred_03 = D2B4_geg[0,3] * (mu(3) - mu(0)) check("P032 [D²,Δ]_geg[0,3] = D²[0,3]·(μ₃-μ₀) = 18432·(-15) = -276480", abs(COM_geg[0,3] - pred_03) < 0.5 and abs(pred_03 - (-276480)) < 0.5) # [1,3]: D²_geg[1,3] * (μ_3 - μ_1) = 18432 * (-12) = -221184 pred_13 = D2B4_geg[1,3] * (mu(3) - mu(1)) check("P033 [D²,Δ]_geg[1,3] = D²[1,3]·(μ₃-μ₁) = 18432·(-12) = -221184", abs(COM_geg[1,3] - pred_13) < 0.5 and abs(pred_13 - (-221184)) < 0.5) check("P034 Eigenvalue gap μ₃ - μ₁ = -12 (= -15 - (-3))", mu(3) - mu(1) == -12) check("P035 Eigenvalue gap μ₃ - μ₀ = -15", mu(3) - mu(0) == -15) check("P036 Eigenvalue gap μ₂ - μ₀ = -8", mu(2) - mu(0) == -8) check("P037 COM_geg has exactly 3 non-zero entries (sparse structure)", sum(1 for v in COM_geg.flat if abs(v) > 1e-6) == 3) # ───────────────────────────────────────────────────────────────────────────── # P038–P042 Layer descent: [D²,Δ]ρ ∈ span{1,x} # ───────────────────────────────────────────────────────────────────────────── print("\n=== P038–P042 Layer descent: [D²,Δ]ρ ∈ span{1,x} ===") # [D²,Δ]ρ in monomial: constant term = 4608π²(π-1), linear = -221184π³ C0_expected = 4608 * PI**2 * (PI - 1) C1_expected = -221184 * PI**3 check("P038 [D²,Δ]ρ has zero x² coefficient (output in span{1,x})", abs(COM_rho[2]) < 1e-6) check("P039 [D²,Δ]ρ has zero x³ coefficient (output in span{1,x})", abs(COM_rho[3]) < 1e-6) check("P040 [D²,Δ]ρ constant term = 4608π²(π-1) (≈ 97398)", abs(COM_rho[0] - C0_expected) < 1.0) check("P041 [D²,Δ]ρ linear coefficient = -221184π³ (≈ -6.858e6)", abs(COM_rho[1] - C1_expected) < 10.0) # The linear coefficient = -221184 * OMEGA_0 * (1/x coefficient) # -221184π³ = -221184 * 4 * OMEGA_0 = -884736 * OMEGA_0 ... let's just verify the exact value check("P042 Linear coeff / constant coeff = -48π/(π-1) (exact ratio)", abs(COM_rho[1]/COM_rho[0] - (-48*PI/(PI-1))) < 1e-8) # ───────────────────────────────────────────────────────────────────────────── # P043–P046 Robertson bound: ratio far from unity # ───────────────────────────────────────────────────────────────────────────── print("\n=== P043–P046 Robertson bound under Haar measure ===") # Expectations under ρ-weighted Haar measure ⟨A⟩ = (1/Ω) ∫₀¹ Aρ · ρ dx # (Here ⟨A⟩ = ∫(Aρ)ρ/∫ρ for the operator expectation in state ρ) # For Robertson: (ΔA)(ΔB) ≥ ½|⟨[A,B]⟩| # We use simpler: ⟨Af⟩ = (1/Ω)∫₀¹ (Af)(x) ρ(x) dx from numpy.polynomial.legendre import leggauss pts, wts = leggauss(100) xs = 0.5 * pts + 0.5 ws = 0.5 * wts rho_xs = 16*PI**3*xs**3 + 3*PI**2*xs**2 + 2*PI*xs # D²_B⁴ ρ = 576π² + 18432π³ x D2B4_rho_xs = 576*PI**2 + 18432*PI**3 * xs # Δ_S³ ρ = 3π + (-6π + 15π²)x + (-24π² + 168π³)x² + (-240π³)x³ DS3_rho_xs = (3*PI + (-6*PI + 15*PI**2)*xs + (-24*PI**2 + 168*PI**3)*xs**2 + (-240*PI**3)*xs**3) # [D²,Δ] ρ = C0 + C1·x COM_rho_xs = C0_expected + C1_expected * xs D2_exp = np.sum(ws * D2B4_rho_xs * rho_xs) / OMEGA DS3_exp = np.sum(ws * DS3_rho_xs * rho_xs) / OMEGA COM_exp = np.sum(ws * COM_rho_xs * rho_xs) / OMEGA # Robertson ratio |⟨[D²,Δ]⟩| / (2|⟨D²⟩||⟨Δ_S³⟩|) robertson = abs(COM_exp) / (2 * abs(D2_exp) * abs(DS3_exp)) CODATA_ALPHA_INV_F = 137.035999084 DELTA = OMEGA - CODATA_ALPHA_INV_F check("P043 ⟨D²_B⁴⟩_ρ > 0 (positive expectation)", D2_exp > 0) check("P044 ⟨Δ_S³⟩_ρ < 0 (negative expectation, consistent with Laplacian spectrum)", DS3_exp < 0) check("P045 Robertson ratio < 0.1 (far from HUP saturation value 1)", robertson < 0.1) check("P046 |⟨[D²,Δ]⟩|/Δ >> 1 (commutator expectation ~ 10¹⁰ × Δ, no algebraic path to Δ)", abs(COM_exp) / DELTA > 1e9) # ───────────────────────────────────────────────────────────────────────────── # Summary # ───────────────────────────────────────────────────────────────────────────── print(f"\nKey values:") print(f" Ω = {OMEGA:.10f}") print(f" OMEGA_0 = π³/4 = {OMEGA_0:.10f}") print(f" D²_B⁴(x²) → 1: {D2B4_mono[0,2]:.0f} (expect 192)") print(f" D²_B⁴(x³) → x: {D2B4_mono[1,3]:.0f} (expect 1152)") print(f" D²_B⁴(U₃)|_U0: {D2B4_geg[0,3]:.1f} (expect 18432)") print(f" D²_B⁴(U₃)|_U1: {D2B4_geg[1,3]:.1f} (expect 18432)") print(f" [D²,Δ]_geg[0,2]: {COM_geg[0,2]:.1f} (expect -24576)") print(f" [D²,Δ]_geg[0,3]: {COM_geg[0,3]:.1f} (expect -276480)") print(f" [D²,Δ]_geg[1,3]: {COM_geg[1,3]:.1f} (expect -221184)") print(f" [D²,Δ]ρ constant: {COM_rho[0]:.4f} (expect {C0_expected:.4f})") print(f" [D²,Δ]ρ linear coeff: {COM_rho[1]:.4e} (expect {C1_expected:.4e})") print(f" Robertson ratio: {robertson:.6e} (HUP minimum = 1)") print(f" |COM_exp| / Δ: {abs(COM_exp)/DELTA:.3e}") print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL") sys.exit(0 if FAIL == 0 else 1)