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

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

def dft(x):
    M=len(x)
    return [
        sum(x[t]*cmath.exp(-2j*math.pi*k*t/M) for t in range(M))/M
        for k in range(M)
    ]

def c2(x):
    M=len(x)
    return [
        sum(x[t]*x[(t+u)%M] for t in range(M))/M
        for u in range(M)
    ]

def c4(x):
    M=len(x)
    out={}
    for u1 in range(M):
        for u2 in range(M):
            for u3 in range(M):
                out[(u1,u2,u3)] = sum(
                    x[t]*x[(t+u1)%M]*x[(t+u2)%M]*x[(t+u3)%M]
                    for t in range(M)
                )/M
    return out

def k4_tensor(x):
    M=len(x)
    C2=c2(x)
    C4=c4(x)
    K={}
    for u1 in range(M):
        for u2 in range(M):
            for u3 in range(M):
                K[(u1,u2,u3)] = (
                    C4[(u1,u2,u3)]
                    - C2[u1]*C2[(u3-u2)%M]
                    - C2[u2]*C2[(u3-u1)%M]
                    - C2[u3]*C2[(u2-u1)%M]
                )
    return K

def Khat_at(K,M,k1,k2,k3):
    return sum(
        K[(u1,u2,u3)]*
        cmath.exp(-2j*math.pi*(k1*u1+k2*u2+k3*u3)/M)
        for u1 in range(M)
        for u2 in range(M)
        for u3 in range(M)
    )/(M**3)

def check_finite_diagonal():
    # odd M avoids nonzero order-2 frequencies
    M=9
    random.seed(82)
    x=[random.uniform(-1,1) for _ in range(M)]
    mean=sum(x)/M
    x=[v-mean for v in x]
    a=dft(x)
    K=k4_tensor(x)
    vals=[]
    D=0.0
    for k in range(1,M):
        kh=Khat_at(K,M,k,k,(-k)%M)
        target=-abs(a[k])**4
        vals.append((k,kh,target))
        assert abs(kh-target)<1e-9,(k,kh,target)
        D += -kh.real
    targetD=sum(abs(a[k])**4 for k in range(1,M))
    assert abs(D-targetD)<1e-9
    return D,targetD

def check_D4_C2():
    M=11
    x=[math.cos(2*math.pi*2*t/M)+0.4*math.cos(2*math.pi*3*t/M+0.3) for t in range(M)]
    a=dft(x)
    C=c2(x)
    lhs=sum(abs(a[k])**4 for k in range(1,M))
    rhs=sum(abs(v)**2 for v in C)/M
    assert abs(lhs-rhs)<1e-10,(lhs,rhs)
    return lhs,rhs

def check_single_mode_moments():
    M=100000
    A=1.3
    vals=[2*A*math.cos(2*math.pi*t/M) for t in range(M)]
    for m in [1,2,3,4]:
        empirical=sum(v**(2*m) for v in vals)/M
        target=math.comb(2*m,m)*(A**(2*m))
        assert abs(empirical-target)/target<1e-10,(m,empirical,target)
    return True

def check_degree4_map():
    for kappa in [0.1,0.3,0.8]:
        d=kappa/2
        threshold=4*d
        assert abs(threshold-2*kappa)<1e-14
        s=threshold+0.02
        kp=s/2
        assert kp>kappa
    return True

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("finite_diagonal",check_finite_diagonal())
    print("D4_equals_C2",check_D4_C2())
    print("single_mode_moments",check_single_mode_moments())
    print("degree4_map",check_degree4_map())
    print("source_delimiters",check_source())
    print("PASS")