from pathlib import Path
import re, math, random

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

def omega(H,r):
    return max(1-abs(r)/H,0.0)

def primes_upto(n):
    sieve=[True]*(n+1)
    sieve[0]=sieve[1]=False
    for p in range(2,int(n**0.5)+1):
        if sieve[p]:
            for k in range(p*p,n+1,p):
                sieve[k]=False
    return [i for i,v in enumerate(sieve) if v]

def von_mangoldt_array(n):
    arr=[0.0]*(n+1)
    ps=primes_upto(n)
    for p in ps:
        x=p
        while x<=n:
            arr[x]=math.log(p)
            x*=p
    return arr

def W(x):
    if x<=1 or x>=2:
        return 0.0
    # compact polynomial bump for finite sanity only
    t=(x-1)*(2-x)
    return t*t

def construct_a(N,H):
    maxn=2*N+H+10
    L=von_mangoldt_array(maxn)
    a=[0.0]*(2*N+1)
    for n in range(N+1,2*N):
        s=0.0
        for r in range(-H+1,H):
            m=n+r
            if m>=1:
                s+=omega(H,r)*L[m]
        a[n]=W(n/N)*s
    return a

def check_local_distribution():
    N=2500
    H=180
    a=construct_a(N,H)
    A=sum(a)
    vals=[]
    for d in [1,2,3,5,7,11,17,31,61,101,181]:
        Ad=sum(a[n] for n in range(len(a)) if n%d==0)
        err=abs(Ad-A/d)
        vals.append((d,err,err/N))
        # finite rough sanity: discrepancy is O(N), not growing like A.
        assert err < 3*N,(d,err,A)
    return vals

def check_kernel_sum():
    for H in [10,50,101]:
        s=sum(omega(H,r) for r in range(-H,H+1))
        assert abs(s-H)<1e-12,(H,s)
    return True

def check_exponent_budget():
    for tau in [0.01,0.05,0.1,0.2]:
        eps=0.01
        eta=1/3-tau-eps
        if tau<1/3-eps:
            assert eta>0
    return True

def gamma_coeff(n,C):
    def mobius(k):
        x=k
        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
    return sum(mobius(d) for d in range(1,min(C,n)+1) if n%d==0)

def check_gamma():
    assert gamma_coeff(30,1)==1
    # gamma(n,C) is integer and divisor-defined.
    for n in range(2,50):
        for C in [1,2,5,10]:
            assert isinstance(gamma_coeff(n,C),int)
    return True

def check_sieve_constant():
    # g(p)=1/p => each finite Euler factor exactly 1
    for p in [2,3,5,7,11,101]:
        fac=(1-1/p)/(1-1/p)
        assert abs(fac-1)<1e-15
    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("kernel_sum",check_kernel_sum())
    print("local_distribution",check_local_distribution())
    print("exponent_budget",check_exponent_budget())
    print("gamma",check_gamma())
    print("sieve_constant",check_sieve_constant())
    print("source_delimiters",check_source())
    print("PASS")