#!/usr/bin/env python3
"""
Synthetic finite-packet crosscheck for CSM_RH Paper 02.

This verifies algebraic identities on an artificial Hilbert-valued exponential
packet. It is NOT numerical evidence for RH and does not use zeta-zero data.

No third-party packages are required.
"""
import cmath, csv, math

rhos = [
    0.76 + 14.10j,
    0.62 + 18.20j,
    0.24 + 14.10j,
]
V = [
    [1.0 + 0.1j, 0.4 - 0.3j, -0.2 + 0.6j],
    [-0.8 + 0.2j, -0.25 + 0.15j, 0.18 - 0.45j],
    [0.35 - 0.1j, 0.2 + 0.05j, 0.15 + 0.1j],
]

def inner(a, b):
    return sum(x * y.conjugate() for x, y in zip(a, b))

def norm2(a):
    return inner(a, a).real

def vec_sum(terms):
    out = [0j] * len(V[0])
    for c, vec in terms:
        for k, z in enumerate(vec):
            out[k] += c * z
    return out

G = [[inner(V[i], V[j]) for j in range(len(V))] for i in range(len(V))]
beta_max = max(r.real for r in rhos)
s_star = 1 + 2 * beta_max
max_idx = max(range(len(rhos)), key=lambda i: rhos[i].real)
expected_residue = norm2(V[max_idx])

def F(T):
    return vec_sum([(cmath.exp(r*T), V[i]) for i, r in enumerate(rhos)])

def E_direct(T):
    return math.exp(T) * norm2(F(T))

def E_pair(T):
    total = 0j
    for i, ri in enumerate(rhos):
        for j, rj in enumerate(rhos):
            total += G[i][j] * cmath.exp((1 + ri + rj.conjugate()) * T)
    return total.real

def L_rational(s):
    total = 0j
    for i, ri in enumerate(rhos):
        for j, rj in enumerate(rhos):
            total += G[i][j] / (s - (1 + ri + rj.conjugate()))
    return total

def stable_laplace_integrand(T, s):
    # Remove the maximal growth exp(beta_max*T) before squaring.
    scaled = vec_sum([
        (cmath.exp((r - beta_max) * T), V[i])
        for i, r in enumerate(rhos)
    ])
    decay = cmath.exp(-(s - s_star) * T)
    return decay * norm2(scaled)

def simpson_complex(func, a, b, n):
    if n % 2:
        n += 1
    h = (b-a)/n
    acc = func(a) + func(b)
    for k in range(1, n):
        acc += (4 if k % 2 else 2) * func(a + k*h)
    return acc * h / 3

def L_numeric(s):
    # s_test is chosen 1.5 to the right of s_star; T=30 gives tail << 1e-19
    return simpson_complex(lambda t: stable_laplace_integrand(t, s), 0.0, 30.0, 30000)

rows = []
for T in [0.0, 0.1, 0.25, 0.5, 0.8]:
    d = E_direct(T)
    p = E_pair(T)
    rows.append(("energy_pair_identity", T, d, p, abs(d-p)))

s_test = s_star + 1.5
lr = L_rational(s_test)
ln = L_numeric(s_test)
rows.append(("laplace_real", s_test, lr.real, ln.real, abs(lr.real-ln.real)))
rows.append(("laplace_imag", s_test, lr.imag, ln.imag, abs(lr.imag-ln.imag)))

# Residue estimate at the maximal diagonal pole.
eps = 1e-8
res_est = eps * L_rational(s_star + eps)
rows.append(("max_pole_residue_real", s_star, expected_residue, res_est.real, abs(expected_residue-res_est.real)))

with open("hlpg_finite_packet_crosscheck.csv", "w", encoding="utf-8", newline="") as f:
    w = csv.writer(f)
    w.writerow(["test", "parameter", "expected_or_direct", "computed", "abs_error"])
    w.writerows(rows)

max_err = max(row[-1] for row in rows)
print(f"beta_max = {beta_max}")
print(f"s_star = {s_star}")
print(f"expected positive diagonal residue = {expected_residue}")
print(f"estimated residue = {res_est}")
print(f"max_abs_error = {max_err:.12e}")
print("status = PASS" if max_err < 1e-5 else "status = FAIL")
raise SystemExit(0 if max_err < 1e-5 else 1)