from pathlib import Path
import re, math, cmath, random

HERE=Path(__file__).resolve().parent
PAPER=HERE/"CSM_RH_Paper_73_Central_Energy_Zero_Strip_Equivalence_and_Log_Scale_Zero_Gram_Orthogonality_v0.1_2026-09-09.md"

def W(u):
    # smooth-ish polynomial bump for numeric sanity, not used as proof object
    if u<=1 or u>=2:
        return 0.0
    x=(u-1)*(2-u)
    return x*x

def G(rho,y,M=5000):
    s=0j
    for j in range(M):
        u=1+(j+0.5)/M
        s += W(u)*(u**(rho-1))*cmath.exp(2j*math.pi*y*u)/M
    return s

def check_G_nonzero():
    rho=0.73+14.1j
    vals=[abs(G(rho,y,2500)) for y in [-0.2,-0.05,0,0.07,0.2]]
    assert max(vals)>1e-6
    return vals

def check_finite_gram_orthogonality():
    # Hilbert space approximated by C^3.
    gammas=[1.3,4.7,9.1]
    V=[
        [1+0.2j,0.4-0.1j,0.1j],
        [0.2+0.3j,0.8+0.1j,-0.2j],
        [0.4,0.1+0.2j,0.6]
    ]
    diag=sum(sum(abs(z)**2 for z in v) for v in V)
    T=20000.0
    M=200000
    acc=0.0
    for j in range(M):
        t=(j+0.5)*T/M
        vec=[0j,0j,0j]
        for g,v in zip(gammas,V):
            phase=cmath.exp(1j*g*t)
            for k in range(3):
                vec[k]+=phase*v[k]
        acc+=sum(abs(z)**2 for z in vec)
    mean=acc/M
    assert abs(mean-diag)/diag<2e-3,(mean,diag)
    return mean,diag

def check_exponent_identity():
    for beta in [0.51,0.63,0.8,0.95]:
        s=2*(1-beta)
        assert abs((1-s/2)-beta)<1e-14
    return True

def check_energy_to_scalar():
    # Algebra only: C=H^2/N ||S||^2 <= N H^2 N^-s => ||S||^2<=N^(2-s)
    for s in [0.1,0.4,0.9]:
        lhs_exp=2-s
        scalar_exp=1-s/2
        assert abs(lhs_exp/2-scalar_exp)<1e-14
    return True

def check_termwise_mellin_kernel():
    # Verify ∫ K(n/N) N^{-z-1} dN = n^{-z} Khat(z)
    # numerically for a simple compact polynomial K.
    def K(u):
        return W(u)*(1+0.2*u)
    z=1.4+0.7j
    n=37.0
    M=20000
    # N ranges n/2 to n; midpoint
    a=n/2; b=n
    integ=0j
    for j in range(M):
        N=a+(j+0.5)*(b-a)/M
        integ+=K(n/N)*(N**(-z-1))*(b-a)/M
    khat=0j
    for j in range(M):
        u=1+(j+0.5)/M
        khat+=K(u)*(u**(z-1))/M
    rhs=(n**(-z))*khat
    assert abs(integ-rhs)<2e-8,(integ,rhs)
    return abs(integ-rhs)

def check_source():
    s=PAPER.read_text(encoding="utf-8")
    for pat in [r"(?<!\\)\\\(",r"(?<!\\)\\\)",r"(?<!\\)\\\[",r"(?<!\\)\\\]"]:
        assert re.search(pat,s) is None
    assert s.count("$$")%2==0
    tmp=re.sub(r"\$\$.*?\$\$","",s,flags=re.S)
    assert len(re.findall(r"(?<!\\)\$",tmp))%2==0
    return True

if __name__=="__main__":
    print("G_nonzero",check_G_nonzero())
    print("finite_gram",check_finite_gram_orthogonality())
    print("exponent_identity",check_exponent_identity())
    print("energy_to_scalar",check_energy_to_scalar())
    print("termwise_mellin_error",check_termwise_mellin_kernel())
    print("source_delimiters",check_source())
    print("PASS")