"""n49 — two GR-divergent signatures made precise, per n49_preregistration.md.
A: short-distance dressing table (2D L=96). B: same-parity interaction sign
(1D ring, info-metric cost, jellium least-squares)."""
import numpy as np

MU2 = 1e-6

# ---------- Part A: dressed vs bare potential, 2D L=96 --------------------
L = 96
k1 = 2*np.pi*np.fft.fftfreq(L)
gr = np.meshgrid(k1, k1, indexing='ij')
kvec = np.stack(gr, axis=-1).reshape(-1, 2)
w = np.sqrt(np.sum(2 - 2*np.cos(kvec.reshape(L, L, 2)), axis=-1).reshape(-1) + MU2)
NS = L*L
coords = np.array(np.unravel_index(np.arange(NS), (L, L))).T
Wq = np.zeros(NS)
for iq in range(NS):
    qc = coords[iq]
    ikp = np.ravel_multi_index(np.mod(coords + qc[None,:], L).T, (L, L))
    kern = 1.0/((w+w[ikp])**2 * w * w[ikp])
    u1 = 1 - np.exp(1j*kvec); u2 = 1 - np.exp(1j*kvec[ikp])
    X = np.conj(u1)*u2
    M = np.einsum('n,na,nb->ab', kern, X.conj(), X)/(8*NS)
    qp = 2*np.pi*qc/L
    Av = 1 + np.exp(-1j*qp)
    try:
        Wq[iq] = np.real(np.conj(Av) @ np.linalg.solve(M, Av))
    except np.linalg.LinAlgError:
        Wq[iq] = np.inf
Wq = Wq.reshape(L, L)
invW = np.where(np.isfinite(Wq) & (np.abs(Wq) > 1e-9), 1/np.where(np.abs(Wq) > 1e-9, Wq, 1), 0.0)
kerW = np.fft.ifftn(invW).real
# bare lattice-Coulomb kernel in the staggered gauge: signless Laplacian
lam_b = 4 + 2*np.cos(gr[0]) + 2*np.cos(gr[1])
inv_b = np.where(lam_b > 1e-9, 1/np.where(lam_b > 1e-9, lam_b, 1), 0.0)
ker_b = np.fft.ifftn(inv_b).real
# potentials: V(r) ~ ker(0)-ker(r); normalize dressed to bare over r in 8..12
rs = np.arange(1, 13)
vd = np.array([kerW[0,0] - kerW[r, 0] for r in rs])
vb = np.array([ker_b[0,0] - ker_b[r, 0] for r in rs])
sc = np.mean(vd[7:12]/vb[7:12])
print("n49-A: short-distance dressing table (2D, on-axis):")
prev = None; mono = True
for i, r in enumerate(rs):
    dev = vd[i]/(sc*vb[i]) - 1
    if r >= 4 and prev is not None and abs(dev) > abs(prev) + 1e-12: mono = False
    if r >= 4: prev = dev
    print(f"  r={r:2d}: deviation from IR Coulomb = {100*dev:+.3f}%")
print(f"  consistency gate (|dev| shrinking for r>=4): {'PASS' if mono else 'FAIL'}")

# ---------- Part B: same-parity pair, 1D info-metric drainage --------------
N = 200

def K_of(dk):
    k = 1.0 + dk
    K = MU2*np.eye(N)
    for e in range(N):
        i, j = e, (e+1) % N
        K[i,i] += k[e]; K[j,j] += k[e]; K[i,j] -= k[e]; K[j,i] -= k[e]
    return K

w2, V = np.linalg.eigh(K_of(np.zeros(N)))
wv = np.sqrt(np.maximum(w2, 1e-30))
U = np.zeros((N, N))
for e in range(N): U[e] = V[e] - V[(e+1) % N]
W1, W2 = np.meshgrid(wv, wv, indexing='ij')
kern = 1.0/((W1+W2)**2 * W1*W2)
B = np.einsum('en,em->enm', U, U).reshape(N, N*N)
G = 0.125*(B*kern.reshape(1, N*N)) @ B.T
A = np.zeros((N, N))
for i in range(N): A[i,(i-1)%N] = 1; A[i,i] += 1
# jellium least-squares drainage (n40 G1 treatment): KKT lstsq with cost G
M = np.zeros((2*N, 2*N)); M[:N,:N] = G; M[:N,N:] = A.T; M[N:,:N] = A
def V_pair(d, F=0.3):
    c = np.zeros(N); c[40] = F; c[(40+d) % N] = F
    rhs = np.zeros(2*N); rhs[N:] = -c
    dk = np.linalg.lstsq(M, rhs, rcond=None)[0][:N]
    return 0.5*dk @ G @ dk
ds = [10,20,30,40,50,60,70,80]
Vs = np.array([V_pair(d) for d in ds])
g_ref = np.array([d*(N-d)/(2*N) for d in ds])
Ad = np.vstack([np.ones_like(g_ref), g_ref]).T
cf, *_ = np.linalg.lstsq(Ad, Vs, rcond=None)
fr = np.sqrt(np.mean((Vs-Ad@cf)**2))/np.ptp(Vs)
print(f"n49-B: same-parity pair: V(d) = {np.round(Vs, 8).tolist()}")
print(f"  slope = {cf[1]:+.5e} -> {'REPULSIVE (as registered: like staggered charges)' if cf[1] < 0 else 'ATTRACTIVE (contradicts Poisson reading)'}")
print(f"  ring-Green fracRMS = {fr:.3%}")