from pathlib import Path
import re, math, cmath

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

def mobius(n):
    x=n
    cnt=0
    p=2
    while p*p<=x:
        if x%p==0:
            x//=p
            cnt+=1
            if x%p==0:
                return 0
        p+=1
    if x>1:
        cnt+=1
    return -1 if cnt%2 else 1

def primes_upto(n):
    out=[]
    for m in range(2,n+1):
        if all(m%d for d in range(2,int(math.sqrt(m))+1)):
            out.append(m)
    return out

def check_centering():
    X,H=43,11
    ps=primes_upto(X)
    C=[sum(mobius(p+h) for p in ps) for h in range(1,H+1)]
    mean=sum(C)/H
    lhs=sum(c*c for c in C)
    rhs=H*mean*mean+sum((c-mean)**2 for c in C)
    assert abs(lhs-rhs)<1e-12
    return lhs,rhs

def check_centered_gram():
    X,H=31,9
    ps=primes_upto(X)
    C=[sum(mobius(p+h) for p in ps) for h in range(1,H+1)]
    mean=sum(C)/H
    V=sum((c-mean)**2 for c in C)
    S={p:sum(mobius(p+h) for h in range(1,H+1)) for p in ps}
    rhs=0.0
    for p1 in ps:
        for p2 in ps:
            k=sum(mobius(p1+h)*mobius(p2+h) for h in range(1,H+1))-S[p1]*S[p2]/H
            rhs+=k
    assert abs(V-rhs)<1e-10
    return V,rhs

def check_fourier_discrete():
    # Cyclic DFT check with L > X+H; the correlation agrees for h<=H.
    X,H=19,6
    L=64
    ps=primes_upto(X)
    a=[0j]*L
    b=[0j]*L
    for p in ps:
        a[p]=1
    for n in range(1,X+H+1):
        b[n]=mobius(n)
    A=[]
    B=[]
    for k in range(L):
        A.append(sum(a[n]*cmath.exp(2j*math.pi*k*n/L) for n in range(L)))
        B.append(sum(b[n]*cmath.exp(2j*math.pi*k*n/L) for n in range(L)))
    errs=[]
    for h in range(1,H+1):
        corr=sum(mobius(p+h) for p in ps)
        four=sum(A[k]*B[k].conjugate()*cmath.exp(2j*math.pi*k*h/L) for k in range(L))/L
        errs.append(abs(corr-four))
        assert abs(corr-four)<1e-10
    return max(errs)

def check_centered_kernel_discrete():
    X,H=17,5
    L=64
    ps=primes_upto(X)
    a=[0j]*L; b=[0j]*L
    for p in ps: a[p]=1
    for n in range(1,X+H+1): b[n]=mobius(n)
    A=[]; B=[]
    for k in range(L):
        A.append(sum(a[n]*cmath.exp(2j*math.pi*k*n/L) for n in range(L)))
        B.append(sum(b[n]*cmath.exp(2j*math.pi*k*n/L) for n in range(L)))
    Z=[A[k]*B[k].conjugate() for k in range(L)]
    C=[sum(mobius(p+h) for p in ps) for h in range(1,H+1)]
    mean=sum(C)/H
    V=sum((c-mean)**2 for c in C)
    # Evaluate Fourier double sum corresponding to centered Dirichlet kernel.
    rhs=0j
    for k in range(L):
        alpha=k/L
        DH_a=sum(cmath.exp(2j*math.pi*h*alpha) for h in range(1,H+1))
        for l in range(L):
            beta=l/L
            DH_diff=sum(cmath.exp(2j*math.pi*h*(alpha-beta)) for h in range(1,H+1))
            DH_b=sum(cmath.exp(2j*math.pi*h*beta) for h in range(1,H+1))
            kernel=DH_diff-DH_a*DH_b.conjugate()/H
            rhs += Z[k]*Z[l].conjugate()*kernel/(L*L)
    assert abs(V-rhs.real)<1e-9 and abs(rhs.imag)<1e-9
    return V,rhs.real

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("centering",check_centering())
    print("centered_gram",check_centered_gram())
    print("fourier_max_error",check_fourier_discrete())
    print("centered_kernel",check_centered_kernel_discrete())
    print("source_delimiters",check_source())
    print("PASS")