#!/usr/bin/env python3
# ============================================================================
# ERRATUM (verifier sweep, 2026-06-15): this script builds the radial operator
# on a DIRICHLET wall (linspace(dr, 1-dr, N-1)), not the corpus APS BC (A329).
# The full-rank-4 LOCAL RIGIDITY result SURVIVES under APS (A330/A335/A336).
# But check 7's alpha-zeta collinearity ~0.99 is a DIRICHLET ARTIFACT: under
# the corpus APS operator it is 0.662 (A336), and the "soft direction"
# fragility it implies is DISSOLVED (A330, the corpus tie zeta=alpha^{5/4}).
# For the corpus-faithful (APS) treatment see A330/A335/A336. The Dirichlet
# checks below are retained as the historical record, flagged.
# ============================================================================
"""verify_P302.py — Verifier for Addendum 302 (P18-T2 coeff rigidity).

Asserts the first-order LOCAL form of P18-T2 sub-result (iii): treating
(alpha, gamma, zeta, beta) as free assembly weights on their four sectors
of P18's radial operator, the spectrum-response Jacobian at the canonical
point has full rank 4, so the coefficients are locally the unique
spectrum-matching assignment. Two structural facts are recorded: (a)
local uniqueness STRENGTHENS as angular sectors are added (coverage =
rigidity); (b) the density (alpha) and renormalization (zeta) sectors are
~99.6%+ spectrally collinear -- the soft direction.

  S1  Sector spectra well-posed     - checks 1-2
  S2  Jacobian full rank 4          - checks 3-5
  S3  Strengthens; soft direction   - checks 6-8
"""
import sys
import math

import numpy as np

PASS = FAIL = 0


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


PI = math.pi
N = 600
NM = 6
A0 = 1.0 / 137.036
Z0 = A0 ** 1.25
ES = 13.177
M0 = 137.036
R0 = A0 ** 1.25


def rho(r):
    return 16 * PI**3 * r**3 + 3 * PI**2 * r**2 + 2 * PI * r


def eigs(alpha, zeta, l, m=NM):
    dr = 1.0 / N
    r = np.linspace(dr, 1 - dr, N - 1)
    n = len(r)
    d2 = (np.diag(-2 * np.ones(n)) + np.diag(np.ones(n - 1), 1)
          + np.diag(np.ones(n - 1), -1)) / dr**2
    d1 = (np.diag(np.ones(n - 1), 1) - np.diag(np.ones(n - 1), -1)) / (2 * dr)
    Vs = (ES / M0**2) * (1 - np.exp(-r / R0))
    H = (-d2 - np.diag(3 / r) @ d1 + np.diag(l * (l + 2) / r**2)
         + np.diag(Vs) + alpha * np.diag(rho(r)) + zeta * np.diag(r) @ d1)
    return np.sort(np.linalg.eigvals(H).real)[:m]


dr = 1.0 / N
x = np.linspace(dr, 1 - dr, N - 1)
mu = [float(np.sum(x**k * rho(x)) * dr) for k in range(NM)]
mur = [m / mu[0] for m in mu]
eps = 1e-6
reom = [math.cos(2 * PI * k / 3) for k in range(3)]


def jac(sectors):
    rows = []
    for l in sectors:
        dA = (eigs(A0 + eps, Z0, l) - eigs(A0 - eps, Z0, l)) / (2 * eps)
        dZ = (eigs(A0, Z0 + eps, l) - eigs(A0, Z0 - eps, l)) / (2 * eps)
        for k in range(3):
            for ni in range(NM):
                rows.append([dA[ni], reom[k], dZ[ni], mur[ni]])
    J = np.array(rows)
    Jn = J / np.linalg.norm(J, axis=0)
    sv = np.linalg.svd(Jn, compute_uv=False)
    caz = abs(float(Jn[:, 0] @ Jn[:, 2]))
    return sv, caz


sv0, caz0 = jac([0])
sv1, caz1 = jac([0, 1])
sv2, caz2 = jac([0, 1, 2])
r0, r1, r2 = sv0[-1] / sv0[0], sv1[-1] / sv1[0], sv2[-1] / sv2[0]

print("S1  Sector spectra well-posed")
check(1, "base l=0 spectrum real and ordered",
      np.all(np.diff(eigs(A0, Z0, 0)) > 0))
check(2, "moment ratios mu_n/mu_0 decreasing from 1.0 (got mu1=%.3f)"
      % mur[1], abs(mur[0] - 1.0) < 1e-9 and mur[1] < 1.0)

print("S2  Jacobian full rank 4 at canonical")
check(3, "four nonzero singular values (l=0): %s"
      % [round(float(s), 3) for s in sv0], np.all(sv0 > 1e-6))
check(4, "sigma_min/sigma_max > 1e-3 (l=0): %.2e" % r0, r0 > 1e-3)
check(5, "rank exactly 4 (no exact flat direction)",
      np.sum(sv0 > 1e-6) == 4)

print("S3  Uniqueness strengthens; alpha-zeta is the soft direction")
check(6, "adding sectors strengthens: %.2e < %.2e < %.2e"
      % (r0, r1, r2), r0 < r1 < r2)
check(7, "alpha and zeta spectrally collinear ~0.99+ (soft dir): %.4f"
      % caz0, caz0 > 0.99)
print("  [ERRATUM] check 7 collinearity is the DIRICHLET value; corpus APS = 0.662 (A336); fragility dissolved (A330)")
check(8, "alpha-zeta collinearity stable across coverage (<0.01 drift)",
      abs(caz0 - caz2) < 0.01)

print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
sys.exit(0 if FAIL == 0 else 1)