# SCRIPT: BH-WALL-03
# Three further checks for the black-hole paper.
#  C9 Where the exceptional points sit. For fixed t and fixed w1-w2 = delta, the shape-block discriminant
#     D(w3) = (A11-A22)^2 + 4 A12 A21 is a quadratic in w3. Claim: every root lies at w3 <= 0
#     (beyond the wall, i.e. imaginary sigma3), and a root reaches w3 = 0 exactly on the joint condition.
#     FALSE if any sampled (t, delta) has a real root in (0, 1-|delta|), or a complex root pair with
#     |Im| > 0 and Re in (0, 1-|delta|) giving D = 0 on the physical segment (checked by direct evaluation).
#  C10 Schur-Cayley boundary. For a random complex L with Herm L >= 0, the Schur complement
#     Z = L_bb - L_bi L_ii^{-1} L_ib satisfies Z = W^* L W with W = [I ; -L_ii^{-1} L_ib], hence
#     Herm Z = W^* Herm L W >= 0, and S = (Z - I)(Z + I)^{-1} has I - S^*S = 4 (Z+I)^{-*} Herm Z (Z+I)^{-1} >= 0.
#     FALSE if any identity is off by > 1e-10 or ||S|| > 1.
#  C11 Local-symbol horizon control. h(r,k) = N(r) (k sigma_x + m sigma_z), N = 1 - 2M/r.
#     Claim: the band projector does not depend on N, so the projector metric in k is the same at every r,
#     including N -> 0; the horizon is invisible to the projector geometry. FALSE if the projector metric at
#     N = 1e-8 differs from N = 1 by > 1e-10.
import numpy as np, sympy as sp
ok = True
t1,t2,t3,w3,dl = sp.symbols('t1 t2 t3 w3 dl', real=True)
w1 = (1 - w3 + dl)/2; w2 = (1 - w3 - dl)/2
A11 = (t1+t2)-(t1-t2)*(w1-w2); A22 = 2*t3 + w3*(t1+t2-2*t3)
A12 = ((t1-t2)-(t1+t2-2*t3)*(w1-w2))/sp.sqrt(3); A21 = sp.sqrt(3)*(t1-t2)*w3
D = sp.expand((A11-A22)**2 + 4*A12*A21)
Dpoly = sp.Poly(D, w3)
print('C9 D(w3) coefficients [w3^2, w3^1, w3^0]:')
for c in Dpoly.all_coeffs(): print('    ', sp.factor(c))
rng = np.random.default_rng(7); worst = -1; bad = 0
f = sp.lambdify((t1,t2,t3,dl), Dpoly.all_coeffs(), 'numpy')
for _ in range(50000):
    t = rng.uniform(0.05, 5, 3); d = rng.uniform(-0.999, 0.999)
    cs = np.array(f(*t, d), dtype=float)
    roots = np.roots(cs) if abs(cs[0]) > 1e-14 else np.roots(cs[1:])
    for r in roots:
        if abs(r.imag) < 1e-12 and 1e-12 < r.real < 1 - abs(d): bad += 1
    if np.all(np.abs(roots.imag) < 1e-12): worst = max(worst, max(roots.real))
print(f'    sampled 50000 (t, delta): roots inside the physical segment: {bad}; largest real root seen: {worst:.3e}')
ok &= bad == 0
c0 = sp.factor(Dpoly.all_coeffs()[-1]); print('    D at the wall (w3=0) =', c0, '  -> zero exactly on the joint condition')
c1 = sp.factor(Dpoly.all_coeffs()[-2])
print('    slope dD/dw3 at w3 = 0 on the joint condition:', sp.factor(sp.simplify(c1.subs(t3, (t1+t2-(t1-t2)*dl)/2))))

rng = np.random.default_rng(3); worstS = 0; err = 0
for _ in range(200):
    n, nb = 6, 2
    X = rng.normal(size=(n,n)) + 1j*rng.normal(size=(n,n)); Hm = X@X.conj().T/n
    K = rng.normal(size=(n,n)) + 1j*rng.normal(size=(n,n)); L = Hm + (K - K.conj().T)/2   # Herm L = Hm >= 0
    Lbb, Lbi, Lib, Lii = L[:nb,:nb], L[:nb,nb:], L[nb:,:nb], L[nb:,nb:]
    Z = Lbb - Lbi@np.linalg.solve(Lii, Lib)
    W = np.vstack([np.eye(nb), -np.linalg.solve(Lii, Lib)])
    err = max(err, np.max(np.abs(Z - W.conj().T@L@W)))
    HZ = (Z + Z.conj().T)/2; err = max(err, np.max(np.abs(HZ - W.conj().T@((L+L.conj().T)/2)@W)))
    S = (Z - np.eye(nb))@np.linalg.inv(Z + np.eye(nb)); Ri = np.linalg.inv(Z + np.eye(nb))
    err = max(err, np.max(np.abs(np.eye(nb) - S.conj().T@S - 4*Ri.conj().T@HZ@Ri)))
    worstS = max(worstS, np.linalg.norm(S, 2))
print(f'C10 Schur-Cayley: max identity error {err:.2e}; max ||S|| over 200 passive samples {worstS:.6f}')
ok &= err < 1e-10 and worstS <= 1 + 1e-12
Z = np.eye(2); S = (Z-np.eye(2))@np.linalg.inv(Z+np.eye(2)); print('     matched load Z = I gives S =', S.ravel())

sx = np.array([[0,1],[1,0]]); sz = np.array([[1,0],[0,-1]])
def proj_metric(N, k, m=1.0, h=1e-5):
    def P(kk):
        H = N*(kk*sx + m*sz); e, V = np.linalg.eigh(H); v = V[:, 1:2]; return v@v.conj().T
    dP = (P(k+h) - P(k-h))/(2*h); return 0.5*np.trace(dP@dP).real
for k in [0.0, 0.7]:
    g1, g0 = proj_metric(1.0, k), proj_metric(1e-8, k)
    print(f'C11 k={k}: projector metric g_kk at N=1: {g1:.12f}; at N=1e-8: {g0:.12f}; exact m^2/(4(k^2+m^2)^2) = {1/(4*(k*k+1)**2):.12f}')
    ok &= abs(g1 - g0) < 1e-10
print('     local proper-time gap 2 sqrt(k^2+m^2) stays open; coordinate gap 2 N sqrt(k^2+m^2) -> 0 is the lapse, not the projector')
print('RESULT', 'PASS' if ok else 'FAIL')