#!/usr/bin/env python3
"""verify_P346b.py -- Verifier for Addendum 346 (P18-T2 sub-result (iii), ANGLE 2 of
a 3-angle attack on the lone non-floor gap: free-weight coefficient->spectrum GLOBAL
injectivity on the corpus APS operator. ANGLE 2 = the OFFSET (+) SHAPE separation:
peel zeta off as a spectral offset, peel gamma off as an exact character block, and
ask whether the remaining (alpha,beta) SHAPE map -- the soft zeta now removed -- is a
P-matrix / globally injective where the full 4x4 (A340/A341/A342) was not.

Copyright Léon Fernando Vlegels -- CC BY 4.0

WHY THIS PROBE EXISTS. A341 (genuine eigh diagonalization of the corpus APS
operator-with-potential H(c)=H0+alpha*rho+zeta*R_sym+beta*M+V_self) found the 4x4
difference-functional Jacobian J_F is NOT a P-matrix (negative principal minor at
every box point) but has a CONSTANT-SIGN, bounded-away-from-zero determinant =>
LOCAL injectivity everywhere; A342 found the classical total-positivity route DEAD
(the kernel is not TP) and controlled the 4x4 det sign only down to one numerical
mixing bound. Both attacked the FULL 4x4 with the soft zeta column IN. ANGLE 2 asks
a structurally different question: zeta is the only non-sign-definite sector (A336);
gamma's character column is exactly c-independent (A342). So SEPARATE them -- recover
zeta from a spectral OFFSET, recover gamma from one character functional, and test
whether the remaining (alpha,beta) SHAPE map closes as a P-matrix on its own. If the
shape block is a clean P-matrix once the soft zeta is removed, that is a route the
full-4x4 P-matrix attack (A340/A341) could not reach.

THE GENUINE OPERATOR (reused VERBATIM from A341/A342 -- the (a) fix). APS basis
R_nl=r^{-1}J_{l+1}(k_n r), k_n=2n+l+2 (all satisfy the APS BC; linear); build the
Galerkin blocks of H(c) in the corpus r^3 measure; mass-orthonormalize via Cholesky;
keep H0=diag(k_n^2) as the corpus APS operator's exact action (A329/P18 l.171); R_sym
from its exact IBP block -2 G + (1/2) psi1 psi1^T; eigh H(c). V=0 HARD GATE:
spectrum (2n+l+2)^2 every l-sector (NOT Dirichlet A333, NOT a shift model A340).

THE FUNCTIONALS (A339's choice, carried). Reference mode (1,2); offset = sum of the
four SEL eigenvalues lambda_(1,2)+lambda_(0,0)+lambda_(2,2)+lambda_(1,1); the three
shape functionals F1=lambda_(0,0)-ref, F2=lambda_(2,2)-ref, F3=lambda_(1,1)-ref.
m=(l,n) -> the n-th genuine (mixed) eigenstate of l-sector l.

WHAT ANGLE 2 ESTABLISHES (the findings, all on the TRUE eigh operator).

  Q1  ZETA OFFSET-RECOVERABLE. d(offset)/d zeta = sum over SEL modes of the TRUE
      <psi_j|R_sym|psi_j> = +5.555 at canonical, and stays SIGN-CONSTANT and bounded
      away from zero ([5.22,5.75]) across the box. Structurally <psi_j|R_sym|psi_j> =
      -2 + b_j with b_j=(1/2)|psi_j(1)|^2/<psi|psi> >= 0 even on the MIXED eigenstate
      (the A333/A334 exact IBP split survives mixing), so the offset's zeta-slope is
      -8 + sum b_j with sum b_j >= 0 -- bounded, nonvanishing, sign-fixed. The offset
      is STRICTLY MONOTONE in zeta at fixed (alpha,gamma,beta) over the box. So GIVEN
      the shape, zeta is GLOBALLY recovered by a strictly-monotone 1D offset map.

  Q2  THE (alpha,beta) SHAPE MAP IS A P-MATRIX (the key new sub-result). gamma's
      character column is EXACTLY c-independent and only F1 carries a nonzero gamma
      difference (1.5); F2 and F3 carry EXACTLY ZERO gamma. So gamma peels off through
      F1, and the (alpha,beta) shape map is the 2x2 [d(F2,F3)/d(alpha,beta)]. Under a
      functional reorder (F3 vs alpha, F2 vs beta) this 2x2 is a P-MATRIX EVERYWHERE on
      the box (both diagonals sign-definite: drho-difference in [19.0,22.7]>0,
      dmom-difference in [0.0259,0.0276]>0; det sign-constant, |det| in [0.103,0.194]).
      Confirmed P-matrix at every grid point + 500 random points. This is EXACTLY the
      sub-problem A340/A341 could not close WITH zeta in the 4x4: with the soft zeta
      REMOVED, the shape block closes as a genuine P-matrix.

  Q3  NON-MIXING -- gamma TRIANGULAR (exact), zeta TRIANGULAR (approximate). gamma
      peels off EXACTLY: F2,F3 have identically zero gamma-derivative (the character
      value carries no eigenstate content, A336). zeta does NOT peel off exactly: the
      TRUE eigenstates shift with zeta, so F2,F3 carry a nonzero zeta-derivative
      (+0.59,-1.34 on the SEL modes) -- the honest coupling. So the offset(+)shape
      structure is block-LOWER-triangular in gamma but only APPROXIMATELY so in zeta.
      The zeta leak into the shape Jacobian is small (as zeta sweeps the full 20x box
      range the (alpha,beta) 2x2 moves by <= 0.039, ~0.2% of its O(10-20) entries), and
      the full 4x4 Phi-Jacobian has a CONSTANT-SIGN, bounded-away-from-zero determinant
      (|det| in [17.5,21.8]) -- LOCAL injectivity of the whole map everywhere -- but the
      exact zeta-coupling means global injectivity of the full 4 does NOT follow purely
      from block injectivity. The corpus-tie family zeta=alpha^{5/4} removes zeta and the
      tied 3x3 (alpha,gamma,beta) shape map is cleanly sign-constant (|det| in
      [0.155,0.290]) -- the cleanly-closed corpus-posed case.

THE HONEST VERDICT -- PARTIAL (with a genuinely new sub-result). ANGLE 2 delivers
something the full-4x4 attacks (A340/A341/A342) could not: the (alpha,beta) SHAPE map,
once the soft zeta and the character gamma are peeled off, is a P-MATRIX everywhere on
the box (sign-definite diagonals + sign-constant det), and zeta is offset-recoverable
by a strictly-monotone 1D map. So the (alpha,gamma,beta) shape block closes rigorously-
modulo-numerical (P-matrix over a dense grid + random, diagonals sign-definite). BUT
the offset(+)shape decomposition is NOT exactly triangular in zeta -- the true
eigenstate shift makes zeta leak (small but nonzero) into the shape functionals -- so
the global injectivity of the FULL free-weight 4 does NOT close from block injectivity
alone; it rests on the small-coupling fact (zeta leak <= 0.039 vs shape-det >= 0.103)
and the A335/A341 constant-sign full determinant, both NUMERICAL. PARTIAL: the SHAPE
sub-block is closed (new), the zeta-offset is monotone-recoverable (new), but the
exact zeta(<->)shape coupling keeps the full-4 closure numerical -- no fourth overclaim.
The corpus tie zeta=alpha^{5/4} is the cleanly-closed corpus-posed-family case.

  S0  Bessel implementation cross-checked against known zeros            1-2
  S1  GENUINE eigh diagonalization: V=0 HARD GATE + live potential       3-5
  S2  Q1: zeta offset-recoverable (d offset/dz sign-fixed, monotone)     6-8
  S3  Q2: the (alpha,beta) SHAPE map IS a P-matrix (the new sub-result)  9-13
  S4  Q3: non-mixing -- gamma exact-triangular, zeta approx (the leak)  14-16
  S5  full-4 status + corpus-tie clean + the honest PARTIAL verdict     17-19

VERDICT: PARTIAL. On the GENUINE eigh-diagonalized corpus APS operator-with-potential
(V=0 gate PASSED, (2n+l+2)^2 every l-sector, NOT Dirichlet/shift), ANGLE 2 separates
zeta (offset) from (alpha,gamma,beta) (shape). NEW SUB-RESULT: with the soft zeta and
the c-independent character gamma peeled off, the remaining (alpha,beta) SHAPE map is a
P-MATRIX EVERYWHERE on the box (sign-definite diagonals, sign-constant det) -- the
sub-problem A340/A341 could not close WITH zeta in the full 4x4. zeta is OFFSET-
RECOVERABLE: d(offset)/d zeta = -8 + sum b_j (b_j>=0 the A333/A334 IBP boundary part,
surviving mixing) is sign-fixed and bounded away from 0, the offset strictly monotone
in zeta. gamma peels off EXACTLY (F2,F3 have zero gamma-derivative). But the offset(+)
shape decomposition is NOT exactly triangular in zeta -- the true eigenstates shift with
zeta so it leaks (small, <=0.039 vs shape-det >=0.103) into the shape functionals -- so
global injectivity of the FULL free-weight 4 does NOT close from block injectivity
alone; it rests on the small coupling + the A335/A341 constant-sign full determinant,
both NUMERICAL. PARTIAL: the shape block closes (new) and zeta is monotone-recoverable
(new), but the exact zeta-shape coupling keeps the full-4 closure numerical. The corpus
tie zeta=alpha^{5/4} is the cleanly-closed corpus-posed-family case. Honesty over
closure -- PARTIAL is reported because the shape block closes but the zeta coupling does
not vanish; no fourth overclaim on this thread.
"""
import sys
import math
import itertools as it

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
A0 = 1.0 / 137.036
Z0 = A0 ** 1.25            # the corpus tie value zeta = alpha^{5/4} (P18 l.348)
G0, B0 = 0.75, 3 * PI / 20  # canonical gamma, beta
ES = 13.177               # E_self (A330 saturating self-lens amplitude)
M0 = 137.036
R0 = A0 ** 1.25
LS = (0, 1, 2)
NM = 10                   # APS Galerkin basis size per l-sector
REOM = [math.cos(2 * PI * k / 3) for k in range(3)]   # Re(omega^k)=1,-1/2,-1/2
J21 = 5.1356223018

BOX = dict(a=(0.6 * A0, 1.6 * A0), g=(0.5 * G0, 1.5 * G0),
           z=(0.2 * Z0, 4.0 * Z0), b=(0.5 * B0, 2.0 * B0))
CANON = (A0, G0, Z0, B0)


# --- numpy-only Bessel J_nu (series) -- the SAME implementation A329/A334/A341 used
def besselj(nu, x):
    x = np.asarray(x, dtype=float)
    s = np.zeros_like(x)
    for m in range(0, 160):
        s = s + ((-1) ** m) * np.exp(
            (2 * m + nu) * np.log(x / 2.0) - math.lgamma(m + 1)
            - math.lgamma(m + nu + 1))
    return s


def jzeros(nu, count):
    xs = np.linspace(1e-4, 45.0, 200000)
    f = besselj(nu, xs)
    idx = np.where(np.diff(np.sign(f)) != 0)[0]
    out = []
    for i in idx[:count]:
        a, b = xs[i], xs[i + 1]
        for _ in range(80):
            mid = 0.5 * (a + b)
            if besselj(nu, np.array([a]))[0] * besselj(nu, np.array([mid]))[0] <= 0:
                b = mid
            else:
                a = mid
        out.append(0.5 * (a + b))
    return np.array(out)


# ===========================================================================
# THE GENUINE GALERKIN OPERATOR-WITH-POTENTIAL (reused from A341/A342 verbatim).
# ===========================================================================
NGR = 12000
RG = np.linspace(1e-6, 1.0, NGR)
DRG = RG[1] - RG[0]
W = RG ** 3 * DRG                          # corpus r^3 measure weight
RHO = 16 * PI ** 3 * RG ** 3 + 3 * PI ** 2 * RG ** 2 + 2 * PI * RG  # P18/P03 density
VSELF = (ES / M0 ** 2) * (1.0 - np.exp(-RG / R0))                   # saturating lens
MOM = RG.copy()                            # beta-sector moment proxy <r> (A335/A336)


def parts(l, nmax=NM):
    """Per-l Galerkin operator blocks in the mass-orthonormalized APS basis (A341)."""
    ks = np.array([2 * n + l + 2 for n in range(nmax)], dtype=float)
    R = np.array([besselj(l + 1, k * RG) / RG for k in ks])          # (nmax, NGR)
    psi1 = np.array([besselj(l + 1, np.array([k]))[0] for k in ks])  # APS boundary vals
    G = (R * W) @ R.T                                                # Gram (r^3 measure)
    Linv = np.linalg.inv(np.linalg.cholesky(G))                     # mass-orthonorm map

    def tr(M):
        return Linv @ M @ Linv.T                                     # block -> ortho basis

    rho_o = tr((R * (RHO * W)) @ R.T)
    vs_o = tr((R * (VSELF * W)) @ R.T)
    mom_o = tr((R * (MOM * W)) @ R.T)
    rs_o = tr(-2.0 * G + 0.5 * np.outer(psi1, psi1))                 # exact-IBP R_sym
    H0 = np.diag(ks ** 2)                                            # APS operator (A329)
    p1_o = Linv @ psi1
    return dict(H0=H0, rho=rho_o, vs=vs_o, mom=mom_o, rs=rs_o, p1=p1_o, ks=ks)


PARTS = {l: parts(l) for l in LS}


def H_of_c(l, c):
    a, g, z, b = c
    B = PARTS[l]
    H = (B['H0'] + a * B['rho'] + z * B['rs'] + b * B['mom'] + B['vs']
         + g * REOM[l % 3] * np.eye(NM))
    return 0.5 * (H + H.T)


def mode_data(l, n, c):
    """eigh the FULL operator-with-potential; read the TRUE diagonals on the n-th
    MIXED eigenvector (the genuine A341 diagonalization)."""
    H = H_of_c(l, c)
    w, V = np.linalg.eigh(H)               # <-- THE GENUINE DIAGONALIZATION
    a_ = V[:, n]
    B = PARTS[l]
    return dict(eig=float(w[n]),
                rho=float(a_ @ B['rho'] @ a_),
                mom=float(a_ @ B['mom'] @ a_),
                rsym=float(a_ @ B['rs'] @ a_),
                b=0.5 * float(a_ @ B['p1']) ** 2)


SEL = [(1, 2), (0, 0), (2, 2), (1, 1)]     # ref, m1, m2, m3 (A339 choice)


def all_minors(M):
    n = M.shape[0]
    out = {}
    for s in range(1, n + 1):
        for idx in it.combinations(range(n), s):
            out[idx] = float(np.linalg.det(M[np.ix_(idx, idx)]))
    return out


# --- box grid (5^4), reused from A335/A341 ---------------------------------
GR = 5
GRID = [(a, g, z, b)
        for a in np.linspace(*BOX['a'], GR)
        for g in np.linspace(*BOX['g'], GR)
        for z in np.linspace(*BOX['z'], GR)
        for b in np.linspace(*BOX['b'], GR)]


def offset(c):
    """The spectral OFFSET functional: sum of the four SEL eigenvalues."""
    return sum(mode_data(l, n, c)['eig'] for (l, n) in SEL)


def shape_row(l, n, c):
    """The (alpha, gamma, beta) Hellmann-Feynman shape response on the TRUE
    eigenstate -- the zeta column DROPPED (the offset/shape separation)."""
    d = mode_data(l, n, c)
    return np.array([d['rho'], REOM[l % 3], d['mom']])


def JF3(c):
    """3x3 shape Jacobian d(F1,F2,F3)/d(alpha,gamma,beta), F_i = lambda_{m_i}-ref."""
    ref = shape_row(*SEL[0], c)
    return np.array([shape_row(*SEL[i + 1], c) - ref for i in range(3)])


def J2_shape_reordered(c):
    """The (alpha,beta) SHAPE 2x2 from F2,F3 (which carry ZERO gamma), reordered so
    F3->alpha-diagonal, F2->beta-diagonal (the P-matrix witness)."""
    ref = shape_row(*SEL[0], c)
    r2 = shape_row(*SEL[2], c) - ref       # F2 = (2,2)-ref
    r3 = shape_row(*SEL[3], c) - ref       # F3 = (1,1)-ref
    # columns (alpha, beta) = indices 0, 2 of the shape row; row0=F3, row1=F2
    return np.array([[r3[0], r3[2]], [r2[0], r2[2]]])


# ===========================================================================
print("S0  Bessel implementation cross-checked against known zeros")

j0_at_0 = besselj(0, np.array([1e-10]))[0]
z1 = jzeros(1, 3)
check(1, "numpy-only Bessel J_nu correct: J_0(0)=%.4f=1, J_1 zeros=%s "
      "(=3.8317,7.0156,10.1735)"
      % (j0_at_0, [round(z, 4) for z in z1]),
      abs(j0_at_0 - 1.0) < 1e-6
      and np.allclose(z1, [3.8317059702, 7.0155866698, 10.1734681351], atol=1e-4))

z2 = jzeros(2, 2)
check(2, "J_2 first zero j_{2,1}=%.4f (=5.1356) -- the same Bessel series "
      "A329/A334/A341 used; the operator is the corpus APS operator" % z2[0],
      abs(z2[0] - J21) < 1e-3)

# ===========================================================================
print("S1  GENUINE eigh diagonalization: V=0 HARD GATE + live potential")

# (3) HARD GATE: at V=0 eigh(H_0) reproduces the integer ladder (2n+l+2)^2 in EVERY
# l-sector (the corpus APS operator, A329; NOT Dirichlet A333, NOT a shift A340).
gate_all = True
gate_str = {}
for l in LS:
    wl, _ = np.linalg.eigh(PARTS[l]['H0'])
    want = [(2 * n + l + 2) ** 2 for n in range(4)]
    gate_str[l] = [int(round(x)) for x in np.sort(wl)[:4]]
    gate_all = gate_all and np.allclose(np.sort(wl)[:4], want)
dir0 = jzeros(1, 4) ** 2
gate_not_dir = not np.allclose(np.sort(np.linalg.eigh(PARTS[0]['H0'])[0])[:3],
                               dir0[:3], rtol=1e-2)
check(3, "HARD GATE: at V=0 eigh(H_0) reproduces (2n+l+2)^2 in EVERY l-sector "
      "(l=0 %s, l=1 %s, l=2 %s), DISTINCT from the Dirichlet squared-Bessel zeros "
      "(%s) -- the GENUINE corpus APS operator, NOT Dirichlet, NOT a shift"
      % (gate_str[0], gate_str[1], gate_str[2],
         [round(float(x), 2) for x in dir0[:3]]),
      gate_all and gate_not_dir)

# (4) THE OPERATOR HAS A LIVE POTENTIAL AND WE eigh IT (the A341 premise carried).
Hc0 = H_of_c(0, CANON)
offdiag = float(np.max(np.abs(Hc0 - np.diag(np.diag(Hc0)))))
w_full, V_full = np.linalg.eigh(Hc0)
second_comp = float(np.sort(np.abs(V_full[:, 0]))[-2])
check(4, "THE OPERATOR HAS A LIVE POTENTIAL AND WE eigh IT: canonical H(c) has "
      "off-diagonal coupling (max off-diag %.4f > 0) and eigh returns MIXED "
      "eigenvectors (second |component| %.4f > 0) -- a TRUE diagonalization (A341)"
      % (offdiag, second_comp), offdiag > 1e-3 and second_comp > 1e-4)

# (5) the A333/A334 exact IBP split SURVIVES the mixing: on the TRUE eigenstate
# <psi_j|R_sym|psi_j> = -2 + b_j with b_j = (1/2)|psi_j(1)|^2/<psi|psi> >= 0.
rsym_sel = np.array([mode_data(l, n, CANON)['rsym'] for (l, n) in SEL])
b_sel = np.array([mode_data(l, n, CANON)['b'] for (l, n) in SEL])
split_ok = np.allclose(rsym_sel, -2.0 + b_sel, atol=1e-6) and np.all(b_sel >= -1e-9)
check(5, "the A333/A334 exact IBP split SURVIVES the genuine mixing: on the TRUE "
      "eigenstate <psi_j|R_sym|psi_j>=-2+b_j, b_j=(1/2)|psi_j(1)|^2/<psi|psi> >= 0 "
      "(SEL: rsym=%s, b=%s, all b>=0, split residual<1e-6) -- the zeta structural "
      "lever holds on the mixed states"
      % ([round(float(x), 3) for x in rsym_sel],
         [round(float(x), 3) for x in b_sel]), split_ok)

# ===========================================================================
print("S2  Q1: zeta offset-recoverable (d offset/dz sign-fixed, monotone)")


def grad_offset_z(c, h=1e-6):
    cp = list(c); cm = list(c); cp[2] += h; cm[2] -= h
    return (offset(tuple(cp)) - offset(tuple(cm))) / (2 * h)


# (6) d(offset)/d zeta at canonical = sum over SEL of the TRUE <psi_j|R_sym|psi_j>
# = -8 + sum b_j (b_j>=0), so it is bounded and sign-fixed by the IBP split.
dz_canon = grad_offset_z(CANON)
identity_val = -8.0 + float(np.sum(b_sel))
check(6, "Q1: d(offset)/d zeta at canonical = %.4f, matching the analytic value "
      "sum_SEL <psi_j|R_sym|psi_j> = -8 + sum b_j = %.4f (b_j>=0 the IBP boundary "
      "part) -- the offset's zeta-slope is bounded and sign-fixed by the IBP split"
      % (dz_canon, identity_val), abs(dz_canon - identity_val) < 1e-2)

# (7) d(offset)/d zeta is SIGN-CONSTANT and bounded away from 0 across the box.
dzs = np.array([grad_offset_z(c) for c in GRID])
dz_sign_const = (len(set(np.sign(dzs))) == 1)
check(7, "Q1: d(offset)/d zeta is SIGN-CONSTANT (all %s) and bounded away from 0 "
      "([%.3f,%.3f]) across the %d^4 box -- zeta is offset-recoverable everywhere"
      % ("positive" if dzs[0] > 0 else "negative", dzs.min(), dzs.max(), GR),
      dz_sign_const and np.min(np.abs(dzs)) > 0.5)

# (8) the OFFSET is STRICTLY MONOTONE in zeta at fixed (alpha,gamma,beta) over the
# box -- so GIVEN the shape, zeta is recovered by a strictly-monotone 1D map.
mono_z = True
for c in GRID[::7]:
    a, g, z0, b = c
    offs = [offset((a, g, zz, b)) for zz in np.linspace(*BOX['z'], 9)]
    dd = np.diff(offs)
    if not (np.all(dd > 0) or np.all(dd < 0)):
        mono_z = False
check(8, "Q1: the OFFSET is STRICTLY MONOTONE in zeta at fixed (alpha,gamma,beta) "
      "across the sampled box (%s) -- GIVEN the shape, zeta is GLOBALLY recovered by "
      "a strictly-monotone 1D offset map (the zeta-recovery result)" % mono_z, mono_z)

# ===========================================================================
print("S3  Q2: the (alpha,beta) SHAPE map IS a P-matrix (the new sub-result)")

# (9) gamma's character column is EXACTLY c-independent and only F1 carries a nonzero
# gamma difference; F2,F3 carry EXACTLY ZERO gamma (the character value, A336).
gvals = [REOM[l % 3] for (l, n) in SEL]            # ref, m1, m2, m3
gdiffs = [gvals[i + 1] - gvals[0] for i in range(3)]   # F1,F2,F3 gamma differences
check(9, "Q2: gamma's character column is EXACTLY c-independent; the gamma "
      "DIFFERENCES are F1=%.3f (nonzero), F2=%.3f, F3=%.3f (EXACTLY ZERO) -- only F1 "
      "carries gamma, so gamma peels off through F1 and the (alpha,beta) shape map is "
      "the F2,F3 2x2 (A336: the character carries no eigenstate content)"
      % (gdiffs[0], gdiffs[1], gdiffs[2]),
      abs(gdiffs[0]) > 1e-9 and abs(gdiffs[1]) < 1e-12 and abs(gdiffs[2]) < 1e-12)

# (10) the reordered (alpha,beta) 2x2 shape map is a P-MATRIX at canonical.
J2c = J2_shape_reordered(CANON)
mm2 = all_minors(J2c)
is_pm_canon = all(v > 0 for v in mm2.values())
check(10, "Q2: the reordered (alpha,beta) SHAPE 2x2 (F3->alpha-diag, F2->beta-diag) "
      "is a P-MATRIX at canonical: diag (%.3f,%.4f) both > 0, det %.5f > 0 (all "
      "principal minors positive) -- the soft zeta REMOVED, the shape block closes"
      % (J2c[0, 0], J2c[1, 1], np.linalg.det(J2c)), is_pm_canon)

# (11) the (alpha,beta) shape 2x2 is a P-MATRIX EVERYWHERE on the box.
pm_all = True
worst_minor = 1e9
dets2 = []
for c in GRID:
    M = J2_shape_reordered(c)
    mm = all_minors(M)
    lo = min(mm.values())
    worst_minor = min(worst_minor, lo)
    dets2.append(np.linalg.det(M))
    if lo <= 0:
        pm_all = False
dets2 = np.array(dets2)
check(11, "Q2: the (alpha,beta) SHAPE 2x2 is a P-MATRIX at EVERY box point (worst "
      "principal minor over the %d^4 box = %.4f > 0; det sign-constant, |det| in "
      "[%.4f,%.4f]) -- the sub-problem A340/A341 could NOT close WITH zeta in the "
      "full 4x4 now closes as a genuine P-matrix with zeta REMOVED"
      % (GR, worst_minor, np.abs(dets2).min(), np.abs(dets2).max()),
      pm_all and worst_minor > 0)

# (12) the P-matrix diagonals are SIGN-DEFINITE over the box (why it is a P-matrix):
# the rho-difference (alpha-diagonal) and mom-difference (beta-diagonal) are fixed-sign.
d0 = np.array([J2_shape_reordered(c)[0, 0] for c in GRID])   # rho-diff, alpha-diag
d1 = np.array([J2_shape_reordered(c)[1, 1] for c in GRID])   # mom-diff, beta-diag
check(12, "Q2: the SHAPE-2x2 diagonals are SIGN-DEFINITE over the box -- alpha-diag "
      "(rho-difference) in [%.2f,%.2f] > 0, beta-diag (mom-difference) in "
      "[%.4f,%.4f] > 0 -- the structural reason the shape block is a P-matrix "
      "(sign-definite difference responses, the A336 shape sectors)"
      % (d0.min(), d0.max(), d1.min(), d1.max()),
      len(set(np.sign(d0))) == 1 and len(set(np.sign(d1))) == 1
      and d0.min() > 0 and d1.min() > 0)

# (13) 500-point random multistart confirms the shape P-matrix (no grid artifact).
rng = np.random.default_rng(0)
pm_rand = True
worst_rand = 1e9
for _ in range(500):
    c = (rng.uniform(*BOX['a']), rng.uniform(*BOX['g']),
         rng.uniform(*BOX['z']), rng.uniform(*BOX['b']))
    lo = min(all_minors(J2_shape_reordered(c)).values())
    worst_rand = min(worst_rand, lo)
    if lo <= 0:
        pm_rand = False
check(13, "Q2: 500-point random multistart confirms the (alpha,beta) SHAPE 2x2 is a "
      "P-MATRIX everywhere (worst principal minor %.4f > 0) -- the new sub-result is "
      "not a grid artifact [NUMERICAL: dense grid + random, the A335-class evidence]"
      % worst_rand, pm_rand and worst_rand > 0)

# ===========================================================================
print("S4  Q3: non-mixing -- gamma exact-triangular, zeta approx (the leak)")


def shape_zeta_deriv(c, h=1e-6):
    """d(F2,F3)/d zeta on the shape functionals -- the zeta leak into the shape."""
    out = []
    for (l, n) in (SEL[2], SEL[3]):
        cp = list(c); cm = list(c); cp[2] += h; cm[2] -= h
        ref_p = mode_data(*SEL[0], tuple(cp))['eig']
        ref_m = mode_data(*SEL[0], tuple(cm))['eig']
        fp = mode_data(l, n, tuple(cp))['eig'] - ref_p
        fm = mode_data(l, n, tuple(cm))['eig'] - ref_m
        out.append((fp - fm) / (2 * h))
    return np.array(out)


# (14) gamma peels off EXACTLY: F2,F3 have IDENTICALLY ZERO gamma-derivative (their
# eigenstates carry the fixed character value, no eigenstate content) -- the gamma
# block is EXACTLY lower-triangular.
def shape_gamma_deriv(c, h=1e-6):
    out = []
    for (l, n) in (SEL[2], SEL[3]):
        cp = list(c); cm = list(c); cp[1] += h; cm[1] -= h
        ref_p = mode_data(*SEL[0], tuple(cp))['eig']
        ref_m = mode_data(*SEL[0], tuple(cm))['eig']
        fp = mode_data(l, n, tuple(cp))['eig'] - ref_p
        fm = mode_data(l, n, tuple(cm))['eig'] - ref_m
        out.append((fp - fm) / (2 * h))
    return np.array(out)
gleak = max(float(np.max(np.abs(shape_gamma_deriv(c)))) for c in GRID[::13])
check(14, "Q3: gamma peels off EXACTLY -- the F2,F3 shape functionals have "
      "identically zero gamma-derivative (max |dF/d gamma| over sampled box = %.2e "
      "~ 0) -- the offset(+)shape structure is EXACTLY lower-triangular in gamma "
      "(the character carries no eigenstate content)" % gleak, gleak < 1e-3)

# (15) zeta does NOT peel off exactly: the TRUE eigenstates shift with zeta, so F2,F3
# carry a NONZERO zeta-derivative (the honest coupling). Quantify it.
sz = shape_zeta_deriv(CANON)
zleak_max = 0.0
for c in GRID[::13]:
    zleak_max = max(zleak_max, float(np.max(np.abs(shape_zeta_deriv(c)))))
check(15, "Q3: zeta does NOT peel off exactly -- the TRUE eigenstates shift with "
      "zeta, so F2,F3 carry a NONZERO zeta-derivative (canonical %s; max over "
      "sampled box %.3f) -- the offset(+)shape decomposition is only APPROXIMATELY "
      "triangular in zeta (the honest coupling, not hidden)"
      % ([round(float(x), 3) for x in sz], zleak_max), zleak_max > 1e-2)

# (16) the zeta leak is SMALL relative to the shape determinant: as zeta sweeps the
# FULL 20x box range, the (alpha,beta) shape 2x2 moves by <= 0.039, ~0.2% of its
# O(10-20) entries, and stays a P-matrix -- the coupling is bounded and weak.
J2ref = J2_shape_reordered(CANON)
shape_move = 0.0
for z in np.linspace(*BOX['z'], 9):
    shape_move = max(shape_move,
                     float(np.max(np.abs(J2_shape_reordered((A0, G0, z, B0)) - J2ref))))
shape_det_min = float(np.abs(dets2).min())
check(16, "Q3: the zeta leak into the SHAPE 2x2 is SMALL: as zeta sweeps the full "
      "20x box range, the (alpha,beta) Jacobian moves by <= %.4f (~%.2f%% of its "
      "O(10-20) entries) and stays a P-matrix; the shape det stays >= %.4f -- the "
      "coupling is bounded and weak relative to the shape block"
      % (shape_move, 100 * shape_move / 20.0, shape_det_min),
      shape_move < 0.2 and shape_det_min > shape_move)

# ===========================================================================
print("S5  full-4 status + corpus-tie clean + the honest PARTIAL verdict")


def Phi_jac(c, h=1e-6):
    """The FULL 4x4 Jacobian d(offset,F1,F2,F3)/d(alpha,gamma,zeta,beta) on the
    genuine eigenvalue differences."""
    def Fvec(cc):
        ref = mode_data(*SEL[0], cc)['eig']
        return np.array([offset(cc),
                         mode_data(*SEL[1], cc)['eig'] - ref,
                         mode_data(*SEL[2], cc)['eig'] - ref,
                         mode_data(*SEL[3], cc)['eig'] - ref])
    J = np.zeros((4, 4))
    for j in range(4):
        cp = list(c); cm = list(c); cp[j] += h; cm[j] -= h
        J[:, j] = (Fvec(tuple(cp)) - Fvec(tuple(cm))) / (2 * h)
    return J


# (17) the FULL 4x4 Phi-Jacobian has a CONSTANT-SIGN, bounded-away-from-zero
# determinant over the box -- LOCAL injectivity of the whole offset(+)shape map
# everywhere (consistent with A341, via these coordinates).
detsF = np.array([float(np.linalg.det(Phi_jac(c))) for c in GRID[::3]])
full_sign_const = (len(set(np.sign(detsF))) == 1)
check(17, "Q3: the FULL 4x4 Phi-Jacobian (offset,F1,F2,F3) vs (alpha,gamma,zeta,beta) "
      "has a CONSTANT-SIGN, bounded-away-from-zero determinant (all %s, |det| in "
      "[%.2f,%.2f]) over the sampled box -- LOCAL injectivity of the whole map "
      "everywhere (consistent with A341 in these coordinates)"
      % ("negative" if detsF[0] < 0 else "positive",
         np.abs(detsF).min(), np.abs(detsF).max()),
      full_sign_const and np.abs(detsF).min() > 1.0)

# (18) the corpus-tie family zeta=alpha^{5/4} is CLEAN: the tied 3x3 (alpha,gamma,beta)
# shape map has a sign-constant determinant over the box (zeta removed entirely).
def JF3_tied(a, g, b):
    z = a ** 1.25

    def sr(l, n):
        d = mode_data(l, n, (a, g, z, b))
        return np.array([d['rho'], REOM[l % 3], d['mom']])
    ref = sr(*SEL[0])
    return np.array([sr(*SEL[i + 1]) - ref for i in range(3)])
dt = []
for a in np.linspace(*BOX['a'], 4):
    for g in np.linspace(*BOX['g'], 4):
        for b in np.linspace(*BOX['b'], 4):
            dt.append(np.linalg.det(JF3_tied(a, g, b)))
dt = np.array(dt)
tie_clean = (len(set(np.sign(dt))) == 1) and np.abs(dt).min() > 1e-3
check(18, "the corpus-tie family zeta=alpha^{5/4} (P18 l.348) is CLEAN: the tied 3x3 "
      "(alpha,gamma,beta) shape map has a SIGN-CONSTANT determinant (|det| in "
      "[%.4f,%.4f]) over the box -- zeta removed entirely, the cleanly-closed "
      "corpus-posed-family case" % (np.abs(dt).min(), np.abs(dt).max()), tie_clean)

# (19) VERDICT: PARTIAL. The (alpha,beta) SHAPE map closes as a P-matrix (new) and
# zeta is offset-recoverable by a strictly-monotone 1D map (new); gamma peels off
# exactly; but the offset(+)shape decomposition is NOT exactly triangular in zeta (the
# eigenstate-shift leak), so global injectivity of the FULL free-weight 4 does NOT
# close from block injectivity alone -- it rests on the small coupling + the A341
# constant-sign full determinant, both NUMERICAL. No fourth overclaim.
shape_block_closes = is_pm_canon and pm_all and pm_rand
zeta_recoverable = dz_sign_const and mono_z
zeta_not_triangular = zleak_max > 1e-2
verdict_partial = (gate_all and offdiag > 1e-3 and split_ok
                   and shape_block_closes and zeta_recoverable
                   and gleak < 1e-3 and zeta_not_triangular
                   and full_sign_const and tie_clean)
check(19, "VERDICT PARTIAL: ANGLE 2 delivers a NEW sub-result the full-4x4 attacks "
      "(A340/A341/A342) could not -- the (alpha,beta) SHAPE map, soft zeta and "
      "character gamma peeled off, is a P-MATRIX everywhere (sign-definite diagonals "
      "+ sign-constant det), and zeta is OFFSET-RECOVERABLE by a strictly-monotone 1D "
      "map; gamma peels off EXACTLY. BUT the offset(+)shape decomposition is NOT "
      "exactly triangular in zeta (the eigenstate-shift leak %.3f), so the FULL "
      "free-weight 4 closure does NOT follow from block injectivity alone -- it rests "
      "on the small coupling + the A341 constant-sign full det, both NUMERICAL. "
      "PARTIAL: shape block closed (new), zeta monotone-recoverable (new), zeta-shape "
      "coupling keeps the full-4 numerical; corpus tie zeta=alpha^{5/4} cleanly closed"
      % zleak_max, verdict_partial)

# ===========================================================================
print(f"\n{'='*70}")
print(f"RESULT: {PASS} PASS / {FAIL} FAIL")
print("VERDICT: PARTIAL (offset(+)shape separation; SHAPE block closes as a P-matrix,")
print("         a NEW sub-result; zeta offset-recoverable; the zeta<->shape coupling")
print("         keeps the FULL-4 closure numerical). ANGLE 2 of the 3-angle attack on")
print("         the lone non-floor gap in P18-T2 (iii) -- free-weight coefficient->")
print("         spectrum GLOBAL injectivity on the corpus APS operator -- separates")
print("         zeta (a spectral OFFSET) from (alpha,gamma,beta) (SHAPE), on the SAME")
print("         GENUINE eigh-diagonalized operator A341/A342 used (V=0 HARD GATE ->")
print("         (2n+l+2)^2 every l-sector, NOT Dirichlet A333, NOT a shift A340; live")
print("         potential eigh'd). Q1 -- ZETA OFFSET-RECOVERABLE: d(offset)/d zeta =")
print("         sum_SEL <psi_j|R_sym|psi_j> = -8 + sum b_j, b_j=(1/2)|psi_j(1)|^2/")
print("         <psi|psi> >= 0 the A333/A334 IBP boundary part SURVIVING the genuine")
print("         mixing, so the offset's zeta-slope is bounded, sign-fixed (+5.555,")
print("         box-range [5.22,5.75]), and the offset is STRICTLY MONOTONE in zeta at")
print("         fixed (alpha,gamma,beta) -- GIVEN the shape, zeta is recovered by a 1D")
print("         monotone map. Q2 -- THE NEW SUB-RESULT: gamma's character column is")
print("         EXACTLY c-independent and only F1 carries a nonzero gamma difference;")
print("         F2,F3 carry EXACTLY ZERO gamma, so peeling gamma off, the remaining")
print("         (alpha,beta) SHAPE 2x2 (F3->alpha-diag, F2->beta-diag) is a P-MATRIX")
print("         EVERYWHERE on the box -- both diagonals sign-definite (rho-difference")
print("         [19.0,22.7]>0, mom-difference [0.0259,0.0276]>0), det sign-constant")
print("         (|det| [0.103,0.194]), P-matrix at every grid point + 500 random. This")
print("         is EXACTLY the sub-problem A340/A341 could NOT close WITH zeta in the")
print("         full 4x4: with the soft zeta REMOVED, the shape block closes as a")
print("         genuine P-matrix. Q3 -- NON-MIXING: gamma peels off EXACTLY (F2,F3 have")
print("         identically zero gamma-derivative, the character carries no eigenstate")
print("         content); zeta does NOT -- the TRUE eigenstates SHIFT with zeta so F2,F3")
print("         carry a nonzero zeta-derivative (+0.59,-1.34 on the SEL modes), the")
print("         honest coupling. The leak is SMALL (as zeta sweeps the full 20x box the")
print("         shape 2x2 moves <= 0.039, ~0.2%, staying a P-matrix; shape det stays")
print("         >= 0.103 > the leak) and the FULL 4x4 Phi-Jacobian has a CONSTANT-SIGN,")
print("         bounded-away-from-zero determinant (|det| [17.5,21.8]) -- LOCAL")
print("         injectivity of the whole map everywhere -- but the exact zeta-coupling")
print("         means GLOBAL injectivity of the full 4 does NOT follow from block")
print("         injectivity alone. CONSEQUENCE: the SHAPE sub-block (alpha,gamma,beta)")
print("         closes rigorously-modulo-numerical (P-matrix over a dense grid + 500")
print("         random; sign-definite diagonals = the A336 shape sectors made a clean")
print("         P-matrix once zeta is out); zeta is offset-recoverable by a monotone 1D")
print("         map; but the eigenstate-shift coupling keeps the FULL free-weight 4")
print("         closure resting on the small-coupling + A341 constant-sign full")
print("         determinant, both NUMERICAL -- so it is a PARTIAL, not a HIT. The")
print("         corpus tie zeta=alpha^{5/4} (P18 l.348) removes zeta entirely and the")
print("         tied 3x3 (alpha,gamma,beta) shape map is cleanly sign-constant -- the")
print("         cleanly-closed corpus-posed-family case. NET FOR (iii): ANGLE 2 SHARPENS")
print("         the floor by closing the (alpha,gamma,beta) SHAPE block as a P-matrix")
print("         (the residual is now ISOLATED to the zeta<->shape coupling, not the")
print("         whole 4x4) and giving zeta a monotone-recovery handle, but does NOT")
print("         close the full free-weight global injectivity -- the exact zeta")
print("         coupling stays numerical. GATE per object: genuine eigh = REUSED")
print("         (A341/A342); V=0 HARD GATE = PASSED; IBP split survives mixing =")
print("         VERIFIED; zeta offset-recovery = SIGN-FIXED + MONOTONE; SHAPE 2x2")
print("         P-matrix = EVERYWHERE (grid + random, the new sub-result); gamma")
print("         triangular = EXACT; zeta triangular = APPROXIMATE (the leak, honest);")
print("         full-4 = constant-sign det (local injectivity, A341); corpus tie =")
print("         CLEAN; closed-form full closure = OPEN (numerical coupling). VERDICT:")
print("         PARTIAL. A294's re-typing of P18-T2 stands; (i)=A328/A329, (ii)=A310,")
print("         (iii) A335 in-box injectivity + A341 constant-sign det unchanged; A336")
print("         per-axis monotonicity unchanged; no canon .tex edited; no published")
print("         number changes. Honesty over closure -- PARTIAL is reported because the")
print("         shape block closes but the zeta coupling does not vanish; no fourth")
print("         overclaim on this thread.")
sys.exit(0 if FAIL == 0 else 1)
