from pathlib import Path
import re, math
import mpmath as mp

HERE=Path(__file__).resolve().parent
PAPER=HERE/"CSM_RH_Paper_80_Cutoff_Flow_Charge_Conservation_and_Empty_Divisor_Hard_Core_v0.1_2026-09-09.md"
mp.mp.dps=70

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 von_mangoldt(n):
    if n<2:return mp.mpf('0')
    x=n; p=2
    while p*p<=x:
        if x%p==0:
            while x%p==0:x//=p
            return mp.log(p) if x==1 else mp.mpf('0')
        p+=1
    return mp.log(n)

def MU_s(U,s):
    return mp.fsum([mobius(d)/(mp.mpf(d)**s) for d in range(1,U+1)])

def LV_s(V,s):
    return mp.fsum([von_mangoldt(e)/(mp.mpf(e)**s) for e in range(1,V+1)])

def constants(U,V):
    MU=MU_s(U,1)
    JU=mp.fsum([mobius(d)*mp.log(d)/d for d in range(1,U+1)])
    LV=LV_s(V,1)
    return MU,JU,LV

def B(U,s):
    return mp.zeta(s)*MU_s(U,s)-1

def E(U,V,s):
    MU,JU,LV=constants(U,V)
    z=mp.zeta(s); zp=mp.diff(mp.zeta,s)
    return zp*(MU-MU_s(U,s))+LV_s(V,s)*(1-z*MU_s(U,s))+(MU*LV+JU+1)*z

def R(U,V,s):
    return E(U,V,s)-mp.zeta(s)

def check_cutoff_flow_zero():
    rho=mp.zetazero(1)
    vals=[]
    for U in [2,3,5,8,13,21]:
        vals.append(complex(B(U,rho)))
    assert max(abs(z+1) for z in vals)<1e-55
    return vals

def check_cutoff_increment_formula():
    s=mp.mpf('1.41')+0.7j
    U1,U2=5,17
    lhs=B(U2,s)-B(U1,s)
    shell=mp.fsum([mobius(d)/(mp.mpf(d)**s) for d in range(U1+1,U2+1)])
    rhs=mp.zeta(s)*shell
    assert abs(lhs-rhs)<mp.mpf('1e-55')
    return float(abs(lhs-rhs))

def check_empty_divisor_zero():
    rho=mp.zetazero(1)
    U=17
    z=mp.zeta(rho)
    empty=z-1
    nonempty=z*(MU_s(U,rho)-1)
    assert abs(empty+1)<1e-55
    assert abs(nonempty)<1e-55
    return complex(empty),complex(nonempty)

def check_R_entire_one_and_zero():
    rho=mp.zetazero(1)
    vals1=[]
    for p in [6,8,10]:
        eps=mp.mpf(10)**(-p)
        vals1.append(complex(R(9,11,1+eps)))
    assert max(abs(z) for z in vals1)<1e6
    valsz=[]
    for p in [5,7,9]:
        eps=mp.mpf(10)**(-p)
        valsz.append(complex(R(9,11,rho+eps)))
    assert max(abs(z) for z in valsz)<1e6
    return vals1[-2:], valsz[-2:]

def divisors(n):
    out=[]
    for d in range(1,int(math.isqrt(n))+1):
        if n%d==0:
            out.append(d)
            if d*d!=n: out.append(n//d)
    return out

def is_squarefree(n):
    p=2
    x=n
    while p*p<=x:
        if x%(p*p)==0:return False
        p+=1
    return True

def bU_int(k,U):
    return sum(mobius(d) for d in divisors(k) if d<=U)

def check_squarefree_duality():
    checks=[]
    for k in range(2,200):
        if not is_squarefree(k): continue
        for U in [1,2,3,5,7,11,17,29]:
            if U>=k: continue
            rhs=-mobius(k)*sum(mobius(q) for q in divisors(k) if q<k/U)
            lhs=bU_int(k,U)
            assert lhs==rhs,(k,U,lhs,rhs)
            checks.append((k,U,lhs))
    return len(checks)

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("cutoff_zero",check_cutoff_flow_zero())
    print("cutoff_increment_error",check_cutoff_increment_formula())
    print("empty_divisor",check_empty_divisor_zero())
    print("R_entire_samples",check_R_entire_one_and_zero())
    print("squarefree_duality_checks",check_squarefree_duality())
    print("source_delimiters",check_source())
    print("PASS")