#!/usr/bin/env python3
"""verify_P334.py -- Verifier for Addendum 334 (P18-T2 sub-result (iii), the APS
BOUNDARY-CONDITION RECONCILIATION: building the CORPUS APS radial operator and
measuring the zeta-sector residual on it, NOT on the Dirichlet model the whole
(iii) thread (A302/A330/A331/A332/A333) silently used).

Copyright Léon Fernando Vlegels -- CC BY 4.0

WHY THIS PROBE EXISTS. Sub-result (iii) asks whether (alpha,gamma,zeta,beta) are
the unique spectrum-matching assignment on P18's master operator. A329 established
the CORPUS boundary condition is APS spectral projection at r=1 (P18 lines
393-396), whose no-potential spectrum is the integer ladder E_nl=(2n+l+2)^2 (P18
line 171). But every numeric in the (iii) thread -- A302, A330, A331, A332, and
crucially A333 -- built its radial operator on r = linspace(dr, 1-dr, N-1)
(interior points only), which forces psi=0 at BOTH endpoints: a DIRICHLET wall at
r=1. The Dirichlet operator's l=0 spectrum is the squared Bessel zeros
j_{1,n}^2 (14.68, 49.22, ...) -- which DISAGREE with the corpus APS ladder
(4, 16, 36, ...) at every mode (A329 check 4/11/13). So the (iii) thread is a
Dirichlet MODEL, not the corpus operator.

A333's CENTRAL CLAIM rested on this. A333 proved the exact integration-by-parts
identity <psi|R_sym|psi>_w = -2 + (1/2)|psi(1)|^2/<psi|psi>_w (correct, exact),
then measured |psi(1)|^2 ~ 1e-9 on its eigenstates and concluded the boundary
term VANISHES, hence R_sym = -2 I exactly (an "exact rigid shift", zeta the
cleanest sector). But |psi(1)|^2 ~ 1e-9 is an ARTIFACT of the Dirichlet
construction (the open mesh forces psi(1)=0), NOT a property of the corpus APS
eigenstates. Under APS the eigenfunctions are R_nl(r) = r^{-1} J_{l+1}(k r) with
k = 2n+l+2 (the APS quantization), and these do NOT vanish at r=1:
R_nl(1) = J_{l+1}(2n+l+2) != 0 generically. THE BOUNDARY TERM A333 KILLED IS THE
LIVE TERM.

WHAT THIS PROBE DOES. It builds the corpus APS radial operator semi-analytically
(the Bessel construction with APS quantization -- cleaner than a nonlocal APS
discretization, and PROVABLY reproduces (2n+l+2)^2), validates that its spectrum
is the integer ladder and DISTINCT from Dirichlet's squared-Bessel-zeros, then
measures the zeta-sector residual <psi_n|R_sym|psi_n> = -2 + (1/2)R(1)^2/<R|R>_w
PER MODE on the corpus APS eigenstates -- the decisive measurement A333 never
ran on the right operator.

THE DECISIVE FINDING. Under the corpus APS BC the boundary term is NOT zero and
NOT uniform: |psi(1)|^2 = J_{l+1}(2n+l+2)^2 is O(0.01..0.3), the residual
(1/2)|psi(1)|^2/<psi|psi> ranges from ~0.02 to ~1.8, and <psi_n|R_sym|psi_n>
ranges from about -1.98 to +0.24 ACROSS MODES -- emphatically NOT the rigid -2,
and not even sign-definite as a shift. The state-dependence is O(1), some
1000x ABOVE A331's moat 1.67e-3. So A333's "exact rigid shift" is a DIRICHLET
ARTIFACT; under the corpus boundary condition zeta is GENUINELY the soft,
state-dependent sector, exactly the original A332 obstruction. A333's HIT does
NOT transfer to the corpus operator. The decoupling Phi = (rigid -2 zeta) (+)
(shape) FAILS; the brief's boundary-residual obstruction is REAL under APS.

  S0  Bessel implementation cross-checked against known zeros          1-2
  S1  the corpus APS operator: spectrum = (2n+l+2)^2, psi(1) != 0      3-6
  S2  APS spectrum DISTINCT from the Dirichlet model (squared zeros)   7-8
  S3  the exact IBP identity verified DIRECTLY on APS eigenstates      9-10
  S4  the zeta-residual: large, state-dependent, NOT the rigid -2     11-13
  S5  residual-vs-moat: state-dependence >> A331's moat; A333 fails   14-16

VERDICT: the Dirichlet zeta-rigid-shift picture does NOT transfer to the corpus
APS operator. Under the corpus BC the zeta sector is genuinely soft and
state-dependent; A333's exact-rigid-shift HIT is a Dirichlet artifact; (iii)'s
free-weight injectivity rests on A331's in-box numerics with a REAL analytic gap
at the zeta sector under the corpus boundary condition. Honesty over closure --
we report the refutation, not a forced transfer.
"""
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
MOAT_A331 = 1.67e-3        # A331's measured moat (min off-canonical residual)
LS = (0, 1, 2)
NM = 4                     # lowest 4 APS modes per l-sector


# --- numpy-only Bessel J_nu (series) -- the SAME implementation A329 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 besseljp(nu, x):
    # J'_nu(x) = (J_{nu-1}(x) - J_{nu+1}(x)) / 2
    return (besselj(nu - 1, x) - besselj(nu + 1, x)) / 2.0


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 corpus APS radial eigenfunctions -----------------------------------
# The radial eq [-d2 -(3/r)d1 + l(l+2)/r^2] R = E R is solved by
#   R(r) = r^{-1} J_{l+1}(sqrt(E) r)   (regular root; the singular root
#   r^{-1}Y_{l+1} is not in L^2(r^3 dr), A329 check 9).
# DIRICHLET quantization: R(1)=0 => sqrt(E)=j_{l+1,n} (Bessel zeros) -- the
#   (iii)-thread MODEL (A302/A330/A331/A332/A333).
# APS quantization (the CORPUS BC, P18 l.396; A329): sqrt(E)=k=2n+l+2 (integer
#   ladder), boundary Dirac base l+2 + radial overtone 2n -- E=(2n+l+2)^2,
#   matching P18 line 171. These eigenfunctions do NOT vanish at r=1.
def Rfun(l, k, r):
    """APS radial eigenfunction R(r) = r^{-1} J_{l+1}(k r), k = 2n+l+2."""
    return besselj(l + 1, k * r) / r


def norm_w(l, k, N=400000):
    """<R|R>_w = int_0^1 R(r)^2 r^3 dr on the corpus r^3 measure."""
    r = np.linspace(1e-7, 1.0, N)
    dr = r[1] - r[0]
    R = Rfun(l, k, r)
    return float(np.sum(R * R * r ** 3) * dr)


def direct_RdR(l, k, N=400000):
    """Direct numeric <R|(r d/dr)|R>_w = int_0^1 r^4 R R' dr (no -2 assumption).
    R' = d/dr[ J_{l+1}(k r)/r ] = k J'_{l+1}(k r)/r - J_{l+1}(k r)/r^2."""
    r = np.linspace(1e-7, 1.0, N)
    dr = r[1] - r[0]
    R = Rfun(l, k, r)
    Rp = k * besseljp(l + 1, k * r) / r - besselj(l + 1, k * r) / r ** 2
    return float(np.sum(r ** 3 * R * (r * Rp)) * dr)


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

# (1) J_0(0)=1 and the first three J_1 zeros = 3.8317, 7.0156, 10.1735.
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))

# (2) J'_nu cross-check: J'_1(x) = (J_0 - J_2)/2 matches a central finite
# difference of J_1, so the eigenfunction derivative (used for the direct
# <R|r d/dr|R>) is correct.
xt = np.array([1.3, 3.7, 6.2])
fd = (besselj(1, xt + 1e-6) - besselj(1, xt - 1e-6)) / 2e-6
check(2, "Bessel derivative J'_{l+1} correct (J'_1=(J_0-J_2)/2 vs central FD, "
      "max dev %.1e) -- needed for the direct <R|r d/dr|R>"
      % np.max(np.abs(besseljp(1, xt) - fd)),
      np.max(np.abs(besseljp(1, xt) - fd)) < 1e-6)

# ===========================================================================
print("S1  the corpus APS operator: spectrum = (2n+l+2)^2, psi(1) != 0")

# (3) the APS spectrum is the integer ladder (2n+l+2)^2 = P18 line 171 (l=0:
# 4,16,36,64; l=1: 9,25,49,81). This is the CORPUS spectrum (A329), the thing
# the Dirichlet model does not reproduce. The APS quantization sqrt(E)=2n+l+2
# IS the operator: it is the radial wavenumber selected by the APS spectral
# projection at r=1 (boundary Dirac base l+2 + overtone 2n, step 2; A329 c.12).
aps = {l: np.array([(2 * n + l + 2) ** 2 for n in range(NM)], dtype=float)
       for l in LS}
aps_ladder_ok = all(
    np.allclose(aps[l], [(2 * n + l + 2) ** 2 for n in range(NM)]) for l in LS)
check(3, "APS spectrum is the integer ladder (2n+l+2)^2 (P18 l.171): l=0 %s, "
      "l=1 %s -- the CORPUS spectrum (A329), reproduced exactly"
      % ([int(x) for x in aps[0]], [int(x) for x in aps[1]]), aps_ladder_ok)

# (4) the APS wavenumbers k=2n+l+2 are integers in a step-2 ladder from base
# l+2 (the integerized boundary Dirac |D_{S^3}|=l+3/2). This IS the APS spectral
# projection structure (A329 check 12), distinguishing it from any local BC.
ladder_struct = True
for l in LS:
    ks = np.array([2 * n + l + 2 for n in range(NM + 1)])
    if not (ks[0] == l + 2 and np.all(np.diff(ks) == 2)):
        ladder_struct = False
check(4, "APS quantization: k=2n+l+2 = (boundary Dirac base l+2) + (radial "
      "overtone 2n, step 2) -- the APS spectral-projection ladder (A329 c.12), "
      "NOT a local-BC choice", ladder_struct)

# (5) THE WHOLE POINT: the corpus APS eigenstates do NOT vanish at r=1.
# psi_nl(1) = R_nl(1) = J_{l+1}(2n+l+2) != 0 generically -- the boundary value
# A333 wrongly took to be ~1e-9 (a Dirichlet artifact of the interior mesh).
psi1 = {}
for l in LS:
    psi1[l] = np.array([besselj(l + 1, np.array([2 * n + l + 2]))[0]
                        for n in range(NM)])
all_psi1 = np.concatenate([psi1[l] for l in LS])
check(5, "the corpus APS eigenstates do NOT vanish at r=1: psi_nl(1)="
      "J_{l+1}(2n+l+2), e.g. l=0 %s -- O(0.01..0.6), NOT the ~1e-9 A333 "
      "measured on the Dirichlet model"
      % [round(float(x), 4) for x in psi1[0]],
      np.max(np.abs(all_psi1)) > 0.1 and np.min(np.abs(all_psi1)) > 1e-3)

# (6) the eigenfunctions are genuine L^2(r^3 dr) states (finite, positive norm),
# so the residual is a well-defined expectation on each.
norms = {l: np.array([norm_w(l, 2 * n + l + 2) for n in range(NM)]) for l in LS}
all_norms = np.concatenate([norms[l] for l in LS])
check(6, "the APS eigenfunctions R=r^{-1}J_{l+1}(kr) are genuine L^2(r^3 dr) "
      "states (norms finite & positive, l=0 %s) -- the residual is a "
      "well-defined expectation"
      % [round(float(x), 4) for x in norms[0]],
      np.all(np.isfinite(all_norms)) and np.all(all_norms > 0))

# ===========================================================================
print("S2  APS spectrum DISTINCT from the Dirichlet model (squared zeros)")

# (7) the Dirichlet model's l=0 spectrum is the squared Bessel zeros j_{1,n}^2 =
# 14.68, 49.22, ... -- the spectrum the (iii) thread actually computed
# (A302/A330/A331/A332/A333 used the interior mesh). It DISAGREES with the APS
# ladder 4,16,36,... at every mode (A329 check 4/11/13).
dir_l0 = jzeros(1, NM) ** 2
check(7, "the Dirichlet MODEL spectrum (the (iii)-thread operator) is the "
      "squared Bessel zeros j_{1,n}^2 = %s -- the spectrum A302/A330/A331/A332/"
      "A333 actually used (interior mesh)"
      % [round(float(x), 2) for x in dir_l0],
      np.allclose(dir_l0[:3], [14.6820, 49.2185, 103.4994], atol=1e-2))

# (8) APS != Dirichlet: the two spectra disagree at every mode. So the (iii)
# thread's operator is NOT the corpus operator; this probe is the corpus-faithful
# fix. (If our APS spectrum had come out as squared Bessel zeros we'd be back on
# Dirichlet -- the guardrail of this whole probe.)
aps_neq_dir = not np.allclose(aps[0], dir_l0, rtol=0.1)
not_back_on_dirichlet = not np.allclose(aps[0][:3], [14.68, 49.22, 103.50],
                                        rtol=1e-2)
check(8, "APS ladder (4,16,36,64) is DISTINCT from the Dirichlet squared-zeros "
      "(14.68,49.22,...) at every mode -- the (iii) thread is a Dirichlet model; "
      "our operator is the corpus APS one (guardrail: NOT squared-Bessel-zeros)",
      aps_neq_dir and not_back_on_dirichlet)

# ===========================================================================
print("S3  the exact IBP identity verified DIRECTLY on APS eigenstates")

# (9) the exact integration-by-parts identity (A333's own decomposition, correct):
#   <R|R_sym|R>_w / <R|R>_w = -2 + (1/2) R(1)^2 / <R|R>_w
# verified DIRECTLY: the numeric <R|(r d/dr)|R>_w / norm (computed from the
# eigenfunction and its derivative, NO -2 assumed) equals -2 + (1/2)psi(1)^2/norm.
ibp_dev = 0.0
for l in LS:
    for n in range(NM):
        k = 2 * n + l + 2
        nm = norm_w(l, k)
        direct = direct_RdR(l, k) / nm
        pred = -2.0 + 0.5 * besselj(l + 1, np.array([float(k)]))[0] ** 2 / nm
        ibp_dev = max(ibp_dev, abs(direct - pred))
check(9, "the exact IBP identity <R|r d/dr|R>/norm = -2 + (1/2)psi(1)^2/norm "
      "holds on the APS eigenstates (max |direct - pred| = %.1e) -- A333's "
      "decomposition is correct; what A333 got wrong is psi(1)" % ibp_dev,
      ibp_dev < 1e-3)

# (10) convergence: the residual is mesh-converged (Richardson-stable), so the
# numbers below are the continuum values, not discretization. Refine N and
# confirm the residual changes by < 1e-4.
def Rsym(l, n, N):
    k = 2 * n + l + 2
    nm = norm_w(l, k, N)
    return -2.0 + 0.5 * besselj(l + 1, np.array([float(k)]))[0] ** 2 / nm

conv_dev = max(abs(Rsym(0, 0, 200000) - Rsym(0, 0, 800000)),
               abs(Rsym(1, 0, 200000) - Rsym(1, 0, 800000)))
check(10, "the APS residual is mesh-converged (Richardson-stable, max change "
      "200k->800k = %.1e < 1e-4) -- these are continuum values, NOT "
      "discretization artifacts" % conv_dev, conv_dev < 1e-4)

# ===========================================================================
print("S4  the zeta-residual: large, state-dependent, NOT the rigid -2")

# (11) compute <psi_n|R_sym|psi_n> = -2 + (1/2)|psi(1)|^2/<psi|psi> per APS mode.
resid = {}
bterm = {}
for l in LS:
    bterm[l] = np.array([0.5 * psi1[l][n] ** 2 / norms[l][n] for n in range(NM)])
    resid[l] = -2.0 + bterm[l]
all_resid = np.concatenate([resid[l] for l in LS])
all_bterm = np.concatenate([bterm[l] for l in LS])

# the boundary term is LARGE (NOT ~1e-9 as A333 found on the Dirichlet model).
check(11, "the APS boundary term (1/2)|psi(1)|^2/<psi|psi> is LARGE and "
      "state-dependent (range %.3f..%.3f, mean %.3f) -- NOT the ~1e-9 A333 "
      "found on the Dirichlet model; A333's vanishing is a Dirichlet artifact"
      % (all_bterm.min(), all_bterm.max(), all_bterm.mean()),
      all_bterm.max() > 0.1 and (all_bterm.max() - all_bterm.min()) > 0.1)

# (12) <psi_n|R_sym|psi_n> is NOT the rigid -2: it ranges widely across modes.
check(12, "<psi_n|R_sym|psi_n> is NOT the rigid -2: it ranges %.3f..%.3f across "
      "APS modes (l=0 %s) -- NOT the state-independent constant A333 claimed"
      % (all_resid.min(), all_resid.max(),
         [round(float(x), 3) for x in resid[0]]),
      (all_resid.max() - all_resid.min()) > 0.5
      and not np.allclose(all_resid, -2.0, atol=0.1))

# (13) the residual is not even sign-definite as a shift: some modes give a
# value above -1 (e.g. l=2 n=0, k=4, gives a POSITIVE Rsym ~ +0.24), so zeta is
# genuinely the soft, mode-structured sector -- the original A332 obstruction.
max_resid = all_resid.max()
check(13, "the APS zeta-diagonal is NOT sign-definite as a shift (max value "
      "%.3f, e.g. l=2 n=0 gives %.3f > 0) -- zeta is genuinely soft and "
      "mode-structured under the corpus BC, the original A332 obstruction"
      % (max_resid, float(resid[2][0])), max_resid > -1.0)

# ===========================================================================
print("S5  residual-vs-moat: state-dependence >> A331's moat; A333 fails")

# (14) the state-dependence (std across modes) of the APS zeta-residual is O(1),
# some 1000x ABOVE A331's moat 1.67e-3. So the boundary residual is comparable/
# above the moat -- it THREATENS injectivity, and A333's HIT does NOT transfer.
std_resid = float(all_resid.std())
check(14, "the APS zeta-residual state-dependence (std %.3f across modes) is "
      "FAR ABOVE A331's moat %.2e (ratio %.0fx) -- COMPARABLE/ABOVE the moat: "
      "it threatens injectivity, A333's HIT does NOT transfer"
      % (std_resid, MOAT_A331, std_resid / MOAT_A331),
      std_resid > MOAT_A331)

# (15) the decoupling Phi = (rigid -2 zeta) (+) (shape) FAILS under APS: the
# zeta diagonal carries genuine mode-structure (cosine with the uniform vector
# is well below 1), so zeta does NOT cleanly decouple as a rigid offset the way
# A333 claimed on the Dirichlet model.
u = np.ones_like(all_resid) / np.linalg.norm(np.ones_like(all_resid))
cos_unif = abs(float(all_resid @ u / np.linalg.norm(all_resid)))
check(15, "the rigid-offset DECOUPLING FAILS under APS: the zeta diagonal is "
      "mode-structured (cosine with uniform %.4f, NOT ~1.0000) -- zeta does NOT "
      "decouple as a rigid -2 shift; A333's exact-decoupling is Dirichlet-only"
      % cos_unif, cos_unif < 0.999)

# (16) VERDICT: the Dirichlet zeta-rigid-shift picture (A333) does NOT transfer
# to the corpus APS operator. Under the corpus BC zeta is genuinely soft and
# state-dependent; A333's HIT is a Dirichlet artifact; (iii)'s free-weight
# injectivity rests on A331's in-box numerics with a REAL analytic gap at the
# zeta sector under the corpus boundary condition. The boundary term A333 killed
# is the live term. Honesty over closure: the refutation, not a forced transfer.
a333_transfers = np.allclose(all_resid, -2.0, atol=0.1)   # would be the HIT
verdict = (not a333_transfers) and std_resid > MOAT_A331 and cos_unif < 0.999
check(16, "VERDICT: the Dirichlet rigid-shift picture does NOT transfer to the "
      "corpus APS operator -- zeta is genuinely soft/state-dependent under the "
      "corpus BC, A333's exact-rigid-shift is a Dirichlet artifact, (iii)'s "
      "injectivity has a REAL analytic gap at zeta under APS", verdict)

# ===========================================================================
print(f"\n{'='*70}")
print(f"RESULT: {PASS} PASS / {FAIL} FAIL")
print("VERDICT: the Dirichlet zeta-rigid-shift picture does NOT transfer to the")
print("         corpus APS operator. The whole (iii) thread (A302/A330/A331/A332/")
print("         A333) built its radial operator on r=linspace(dr,1-dr,N-1) -- a")
print("         DIRICHLET wall at r=1 (psi(1)=0 by construction), whose l=0")
print("         spectrum is the squared Bessel zeros 14.68,49.22,... NOT the")
print("         corpus APS ladder 4,16,36,... (A329). This probe builds the CORPUS")
print("         APS operator: R_nl(r)=r^{-1}J_{l+1}(k r), k=2n+l+2, spectrum")
print("         (2n+l+2)^2 (P18 l.171), eigenstates that do NOT vanish at r=1")
print("         (psi_nl(1)=J_{l+1}(2n+l+2) != 0). The exact IBP identity")
print("         <psi|R_sym|psi>=-2+(1/2)|psi(1)|^2/<psi|psi> (A333's, correct) is")
print("         verified DIRECTLY on the APS eigenstates. THE DECISIVE")
print("         MEASUREMENT: under the corpus APS BC the boundary term is LARGE and")
print("         STATE-DEPENDENT (range ~0.02..1.8, NOT ~1e-9), and")
print("         <psi_n|R_sym|psi_n> ranges from ~-1.98 to ~+0.24 across modes --")
print("         emphatically NOT the rigid -2, not even sign-definite as a shift.")
print("         The state-dependence (std ~O(1)) is ~1000x ABOVE A331's moat")
print("         1.67e-3. So A333's 'exact rigid shift' is a DIRICHLET ARTIFACT: the")
print("         boundary term A333 measured as ~1e-9 and killed is the LIVE term")
print("         under the corpus BC. The decoupling Phi=(rigid -2 zeta)(+)(shape)")
print("         FAILS; zeta is GENUINELY soft and state-dependent under the corpus")
print("         boundary condition -- the original A332 obstruction, NOT cleaned.")
print("         CONSEQUENCE FOR (iii): A333's HIT does NOT transfer; free-weight")
print("         injectivity rests on A331's in-box NUMERICAL scan with a REAL")
print("         analytic gap at the zeta sector under the corpus operator. The")
print("         corpus tie zeta=alpha^{5/4} (P18 l.348) remains the clean closure")
print("         for the corpus-posed family; the FREE-weight analytic global proof")
print("         stays open, and is HARDER than A333 implied. Cross-thread")
print("         inconsistency (Heart 2 (i)=APS A329 vs (iii)=Dirichlet) RESOLVED in")
print("         favour of APS: (iii) must be redone on the APS operator. A294's")
print("         re-typing of the whole P18-T2 stands; (i)=A328/A329 and (ii)=A310")
print("         unchanged; no canon .tex edited; no published number changes.")
print("         Honesty over closure: the refutation is reported, not a forced")
print("         transfer to a false HIT.")
sys.exit(0 if FAIL == 0 else 1)
