from pathlib import Path
import re, math, random

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

def chain_bound(A,H,p):
    N=(len(A)-1)//2
    lhs=sum(abs(A[n])**p for n in range(1,2*N+1))
    anchor=sum(abs(A[r])**p for r in range(1,min(H,2*N)+1))
    lag=0.0
    for n in range(1,2*N-H+1):
        lag += abs(A[n+H]-A[n])**p
    M=math.ceil(2*N/H)+1
    rhs=(2**(p-1))*(M*anchor + (M**p)*lag)
    return lhs,rhs

def check_chain():
    random.seed(59)
    worst=0.0
    for p in [1,1.5,2,3,5]:
        for N in [20,50,100]:
            for H in [2,5,11]:
                if H>=2*N: continue
                A=[0.0]+[random.uniform(-3,3) for _ in range(2*N)]
                lhs,rhs=chain_bound(A,H,p)
                assert lhs<=rhs+1e-9,(p,N,H,lhs,rhs)
                worst=max(worst,lhs/rhs if rhs else 0)
    return worst

def d_p(k,a,nu,c,p):
    d=k/2
    return min(nu,d+a-1+c/p,1-a*(1-d))

def phi1(k,a,nu,c):
    return min(2*nu,k+2*c+2*a-2,2-2*a+a*k)

def check_l1_optimality():
    random.seed(590)
    for _ in range(10000):
        k=random.uniform(0.02,0.98)
        a=random.uniform(0.05,0.99)
        nu=random.uniform(0.01,0.8)
        c=random.uniform(0.01,1.0)
        best=d_p(k,a,nu,c,1)
        for p in [1.1,1.5,2,3,5,10]:
            assert d_p(k,a,nu,c,p)<=best+1e-14
        assert abs(2*best-phi1(k,a,nu,c))<1e-14
    return True

def check_strict_gate():
    random.seed(591)
    for _ in range(20000):
        k=random.uniform(0.02,0.98)
        a=random.uniform(0.02,0.999)
        nu=random.uniform(0.001,0.8)
        c=random.uniform(0.001,1.0)
        amp=phi1(k,a,nu,c)>k+1e-12
        gate=(a<1 and nu>k/2 and c>1-a)
        assert amp==gate,(k,a,nu,c,phi1(k,a,nu,c),amp,gate)
    return True

def check_optimum():
    worst=0.0
    for k in [0.1,0.3,0.5,0.7,0.9]:
        for c in [0.02,0.05,0.1,0.2]:
            a_star=1-2*c/(4-k)
            target=k+2*c*(2-k)/(4-k)
            nu=0.7
            vals=[]
            for j in range(1,40000):
                a=j/40000
                vals.append(phi1(k,a,nu,c))
            num=max(vals)
            # If target exceeds threshold cap or 1, skip exact optimiser comparison.
            expected=min(target,2*nu)
            worst=max(worst,abs(num-expected))
            assert abs(num-expected)<1e-4,(k,c,a_star,num,expected)
    return worst

def check_moment_boundary():
    random.seed(592)
    for _ in range(10000):
        k=random.uniform(0.05,0.95)
        d=k/2
        nu=random.uniform(d+1e-4,min(0.9,d+0.3))
        q=random.uniform(0.01,0.8)
        A=max(0.0,random.uniform(0,3))
        for r in [1,2,3,5,10]:
            exponent=1+2*r*(nu-d)+q*A*d
            assert exponent>1
    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("chain_ratio_worst",check_chain())
    print("L1_optimality",check_l1_optimality())
    print("strict_gate",check_strict_gate())
    print("optimizer_error",check_optimum())
    print("moment_boundary",check_moment_boundary())
    print("source_delimiters",check_source())
    print("PASS")