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

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

def F(rho,x):
    return -(x**rho)/rho

def lag_energy(N,H,rho):
    return sum(abs(F(rho,n+H)-F(rho,n))**2 for n in range(N,2*N-H))

def defect_energy(N,H,rho):
    vals=[F(rho,n+H)-F(rho,n) for n in range(N,2*N-2*H)]
    vals2=[F(rho,n+2*H)-F(rho,n+H) for n in range(N,2*N-2*H)]
    return sum(abs(a-b)**2 for a,b in zip(vals,vals2))

def check_zero_mode_scaling():
    rho=0.7+14.1347251417347j
    alpha=0.65
    ratios=[]
    defects=[]
    for N in [4000,8000,16000,32000]:
        H=max(2,int(N**alpha))
        S=lag_energy(N,H,rho)
        scale=(H**2)*(N**(2*rho.real-1))
        ratios.append(S/scale)
        D=defect_energy(N,H,rho)
        defects.append((D/S)/((H/N)**2))
    assert min(ratios)>0.001 and max(ratios)<1000
    assert min(defects)>0.001 and max(defects)<1e6
    return ratios,defects

def increments(N,rho):
    return [F(rho,n)-F(rho,n-1) for n in range(N+1,2*N+1)]

def fourier_mass_fraction_outside(N,H,rho,c=0.5,M=40000):
    b=increments(N,rho)
    total=sum(abs(z)**2 for z in b)
    # numerical trapezoid/midpoint integration over [0,1), mark distance to nearest integer.
    s=0.0
    dx=1.0/M
    for j in range(M):
        xi=(j+0.5)*dx
        norm=min(xi,1-xi)
        if norm>c/H:
            val=0j
            # use a smaller N for this validation only
            for idx,z in enumerate(b, start=N+1):
                val += z*cmath.exp(2j*math.pi*idx*xi)
            s += abs(val)**2*dx
    return s/total

def check_principal_concentration():
    rho=0.7+14.1347251417347j
    vals=[]
    # modest N to keep validation fast
    for N in [120,180,260]:
        H=max(3,int(N**0.6))
        frac=fourier_mass_fraction_outside(N,H,rho,M=5000)
        vals.append((N,H,frac, H/N))
    # finite-N sanity only: outside fraction should decrease and remain compatible with O(H/N)
    assert vals[-1][2] < vals[0][2]
    assert all(frac < 20*(H/N) for N,H,frac,_ in vals)
    return vals

def check_exponent_law_algebra():
    for k in [0.1,0.25,0.5,0.8,1.0]:
        sigma=1-k/2
        assert abs((2*sigma+1)-(3-k))<1e-14
    return True

def check_source():
    s=PAPER.read_text(encoding="utf-8")
    forbidden=[
        r"(?<!\\)\\\(",
        r"(?<!\\)\\\)",
        r"(?<!\\)\\\[",
        r"(?<!\\)\\\]",
    ]
    for pat in forbidden:
        assert re.search(pat,s) is None, pat
    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("zero_mode_scaling",check_zero_mode_scaling())
    print("principal_concentration",check_principal_concentration())
    print("exponent_law_algebra",check_exponent_law_algebra())
    print("source_delimiters",check_source())
    print("PASS")