#!/usr/bin/env python3
"""Exact scalar interval-positivity and finite-moment certificates.

Standard library only. Gram certificates and rational atoms are supplied,
not found by a general SDP or algebraic-root solver. The terminal quadratic
pair is verified by rational isolating intervals and a power recurrence.
"""
from fractions import Fraction as F
from itertools import product,combinations
from pathlib import Path
from collections import Counter
import argparse,json
COUNTS=Counter()
def check(condition,label):
    COUNTS[label]+=1
    if not condition:raise ArithmeticError(label)
def trim(p):
    p=list(map(F,p))
    while len(p)>1 and not p[-1]:p.pop()
    return p or [F(0)]
def add(p,q):
    return trim([(p[i] if i<len(p) else 0)+(q[i] if i<len(q) else 0) for i in range(max(len(p),len(q)))])
def scale(p,c):return trim([F(c)*x for x in p])
def mul(p,q):
    out=[F(0)]*(len(p)+len(q)-1)
    for i,x in enumerate(p):
        for j,y in enumerate(q):out[i+j]+=x*y
    return trim(out)
def value(p,x):
    y=F(0)
    for c in reversed(p):y=y*x+c
    return y
def L(p,m):
    p=trim(p)
    if len(p)>len(m):raise ValueError('unknown moment needed')
    return sum((x*F(m[i]) for i,x in enumerate(p)),F(0))
def dot(x,y):return sum((a*b for a,b in zip(x,y)),F(0))
def quadratic(M,v):return sum((v[i]*M[i][j]*v[j] for i in range(len(v)) for j in range(len(v))),F(0))
def psd(M):
    """Exact recursive Schur certificate, including zero pivots and a negative vector."""
    M=[list(map(F,row)) for row in M];n=len(M)
    if any(len(row)!=n for row in M) or any(M[i][j]!=M[j][i] for i in range(n) for j in range(n)):
        raise ValueError('real symmetric square matrix required')
    if not n:return {'psd':True,'rank':0,'kernel':[],'negative':None}
    pivot=M[0][0]
    if pivot<0:return {'psd':False,'negative':[F(1)]+[F(0)]*(n-1)}
    b=M[0][1:];D=[row[1:] for row in M[1:]]
    if pivot==0:
        j=next((j for j in range(1,n) if M[0][j]),None)
        if j is not None:
            v=[F(0)]*n;v[j]=1;v[0]=-(M[j][j]+1)/(2*M[0][j])
            return {'psd':False,'negative':v}
        s=psd(D)
        if not s['psd']:return {'psd':False,'negative':[F(0)]+s['negative']}
        return {'psd':True,'rank':s['rank'],'kernel':[[F(1)]+[F(0)]*(n-1)]+[[F(0)]+v for v in s['kernel']],'negative':None}
    S=[[D[i][j]-b[i]*b[j]/pivot for j in range(n-1)] for i in range(n-1)]
    s=psd(S)
    def lift(v):return [-dot(b,v)/pivot]+v
    if not s['psd']:return {'psd':False,'negative':lift(s['negative'])}
    return {'psd':True,'rank':1+s['rank'],'kernel':[lift(v) for v in s['kernel']],'negative':None}
def gram_polynomial(G):
    out=[F(0)]*max(1,2*len(G)-1)
    for i,row in enumerate(G):
        for j,x in enumerate(row):out[i+j]+=F(x)
    return trim(out)
def gram_certificate(p,a,b,n,parity,G0,G1):
    a,b=F(a),F(b)
    if a>=b or not isinstance(n,int) or n<0 or parity not in ('even','odd'):
        raise ValueError('nondegenerate interval, integer n>=0 and even/odd required')
    sizes=(n+1,n if parity=='even' else n+1)
    if len(G0)!=sizes[0] or len(G1)!=sizes[1]:raise ValueError('wrong Gram dimensions')
    if len(trim(p))>2*n+(1 if parity=='even' else 2):raise ValueError('degree exceeds contract')
    if not psd(G0)['psd'] or not psd(G1)['psd']:return False
    if parity=='even':candidate=add(gram_polynomial(G0),mul([-a*b,a+b,-1],gram_polynomial(G1)))
    else:candidate=add(mul([-a,1],gram_polynomial(G0)),mul([b,-1],gram_polynomial(G1)))
    return candidate==trim(p)
def localizers(m,a,b):
    m=list(map(F,m));a,b=F(a),F(b)
    if not m or a>=b:raise ValueError('full nonempty prefix and a<b required')
    d=len(m)-1;n=d//2
    def matrix(weight,size):return [[L([0]*(i+j)+weight,m) for j in range(size)] for i in range(size)]
    if d%2==0:
        return [('H',[F(1)],matrix([F(1)],n+1)),('K',[-a*b,a+b,F(-1)],matrix([-a*b,a+b,F(-1)],n))]
    return [('A',[-a,F(1)],matrix([-a,F(1)],n+1)),('B',[b,F(-1)],matrix([b,F(-1)],n+1))]
def interval_moments(m,a=0,b=1):
    m=list(map(F,m));a,b=F(a),F(b)
    if not m or a>b:raise ValueError('nonempty full prefix and ordered interval required')
    if a==b:
        if m[0]<0:return {'status':'INFEASIBLE','witness':[F(1)]}
        for k,mk in enumerate(m):
            delta=mk-m[0]*a**k
            if delta:
                p=[-a**k]+[F(0)]*(k-1)+[F(1)]
                return {'status':'INFEASIBLE','witness':scale(p,-1 if delta>0 else 1)}
        return {'status':'ZERO' if m[0]==0 else 'UNIQUE','support_polynomial':[-a,F(1)],'rank_records':[]}
    records=[];support=None
    for name,weight,M in localizers(m,a,b):
        cert=psd(M)
        if not cert['psd']:
            q=cert['negative'];witness=mul(weight,mul(q,q))
            if L(witness,m)>=0:raise ArithmeticError('invalid negative certificate')
            return {'status':'INFEASIBLE','matrix':name,'vector':q,'witness':witness}
        records.append((name,len(M),cert['rank']))
        if cert['kernel'] and support is None:support=mul(weight,mul(cert['kernel'][0],cert['kernel'][0]))
    if m[0]==0:
        if any(m):raise ArithmeticError('zero mass inconsistency')
        return {'status':'ZERO','rank_records':records}
    return {'status':'UNIQUE' if support is not None else 'MULTIPLE','support_polynomial':support,'rank_records':records}
def rational_atom_certificate(m,a,b,nodes,weights):
    a,b=F(a),F(b);nodes=list(map(F,nodes));weights=list(map(F,weights))
    if not m or a>b or len(nodes)!=len(weights) or len(set(nodes))!=len(nodes):return False
    if any(x<a or x>b for x in nodes) or any(w<=0 for w in weights):return False
    return all(sum((w*x**k for x,w in zip(nodes,weights)),F(0))==F(mk) for k,mk in enumerate(m))
def atoms(nodes,weights,d):return [sum((F(w)*F(x)**k for x,w in zip(nodes,weights)),F(0)) for k in range(d+1)]
def outer(v):return [[F(x)*F(y) for y in v] for x in v]
def determinant(M):
    # Separate small-matrix Laplace oracle, not the PSD elimination path.
    if not M:return F(1)
    return sum(((-1)**j*M[0][j]*determinant([r[:j]+r[j+1:] for r in M[1:]]) for j in range(len(M))),F(0))
def all_principal_nonnegative(M):
    n=len(M)
    return all(determinant([[M[i][j] for j in inds] for i in inds])>=0 for k in range(1,n+1) for inds in combinations(range(n),k))
def run():
    COUNTS.clear()
    # Exhaustive exact indefinite/PSD/zero-pivot cases with an independent determinant criterion.
    for data in product([-1,0,1],repeat=6):
        x,y,z,u,v,w=map(F,data);M=[[x,y,z],[y,u,v],[z,v,w]];cert=psd(M)
        check(cert['psd']==all_principal_nonnegative(M),'PSD versus all principal minors')
        if cert['psd']:
            check(len(cert['kernel'])==3-cert['rank'],'exact nullity')
            for q in cert['kernel']:check(all(dot(row,q)==0 for row in M),'kernel annihilation')
        else:check(quadratic(M,cert['negative'])<0,'negative quadratic direction')
    # Rational Gram identities at several affine intervals and degree bounds.
    for a,b in [(F(0),F(1)),(F(-1),F(1)),(F(2),F(5)),(F(-3),F(-1,2))]:
        for n in range(4):
            for seed in range(7):
                v=[F((seed+2*j)%5-2,3) for j in range(n+1)]
                G=outer(v)
                for parity in ['even','odd']:
                    sz=n if parity=='even' else n+1;t=[F((seed+3*j)%7-3,4) for j in range(sz)];D=outer(t)
                    if parity=='even':p=add(mul(v,v),mul([-a*b,a+b,-1],mul(t or [0],t or [0])))
                    else:p=add(mul([-a,1],mul(v,v)),mul([b,-1],mul(t,t)))
                    check(gram_certificate(p,a,b,n,parity,G,D),'supplied interval Gram coefficient identity')
                    check(not gram_certificate(add(p,[F(1,97)]),a,b,n,parity,G,D),'one coefficient mutation rejected')
                    for k in range(6):check(value(p,a+(b-a)*F(k,5))>=0,'supplementary rational value checks')
    # Moment matrices from supplied positive atoms; rank/uniqueness uses the theorem.
    for a,b in [(F(0),F(1)),(F(-2),F(3)),(F(2),F(4))]:
        grid=[a+(b-a)*F(k,4) for k in range(5)]
        for subset_size in range(1,5):
            for inds in combinations(range(5),subset_size):
                xs=[grid[i] for i in inds];weights=[F(j+1) for j in range(subset_size)]
                for d in range(9):
                    m=atoms(xs,weights,d);cert=interval_moments(m,a,b)
                    check(cert['status'] in ('UNIQUE','MULTIPLE'),'positive atom prefix feasibility')
                    check(rational_atom_certificate(m,a,b,xs,weights),'supplied rational atoms all moments')
                    for name,weight,M in localizers(m,a,b):
                        expected=[[sum((w*value(weight,x)*x**(i+j) for x,w in zip(xs,weights)),F(0)) for j in range(len(M))] for i in range(len(M))]
                        check(M==expected,'localizer integral identity')
                    if cert['status']=='UNIQUE':check(all(value(cert['support_polynomial'],x)==0 for x in xs),'unique support zeros contain all atoms')
    # Four terminal endpoints and their negative perturbations.
    prefix=[F(1),F(1,2),F(1,3)];third=[]
    for t,xs,ws in [(F(2,9),[0,F(2,3)],[F(1,4),F(3,4)]),(F(5,18),[F(1,3),1],[F(3,4),F(1,4)])]:
        m=prefix+[t];check(interval_moments(m)['status']=='UNIQUE','third-moment endpoint uniqueness');check(rational_atom_certificate(m,0,1,xs,ws),'third-moment endpoint measure');third.append(m)
    for k in range(21):
        t=F(2,9)+F(k,20)*(F(5,18)-F(2,9));check(interval_moments(prefix+[t])['status']==('UNIQUE' if k in (0,20) else 'MULTIPLE'),'complete third-moment interval')
    for t in [F(2,9)-F(1,1000),F(5,18)+F(1,1000)]:check(interval_moments(prefix+[t])['status']=='INFEASIBLE','outside third-moment interval')
    q=[F(1,6),F(-1),F(1)];g=[F(0),F(1),F(-1)];ell=[F(-1,2),F(1)]
    lower=mul(q,q);upper=mul(g,mul(ell,ell));fourth_prefix=prefix+[F(1,4)]
    for k in range(25):
        s=F(7,36)+F(k,24)*(F(5,24)-F(7,36));m=fourth_prefix+[s]
        check(L(lower,m)==s-F(7,36) and L(upper,m)==F(5,24)-s,'fourth-moment dual polynomial identities')
        check(interval_moments(m)['status']==('UNIQUE' if k in (0,24) else 'MULTIPLE'),'complete fourth-moment interval')
    upper_m=fourth_prefix+[F(5,24)];check(rational_atom_certificate(upper_m,0,1,[0,F(1,2),1],[F(1,6),F(2,3),F(1,6)]),'upper endpoint three-node measure')
    loc=localizers(upper_m,0,1);check(psd(loc[0][2])['rank']==3 and psd(loc[1][2])['rank']==1,'positive Hankel with singular localizer')
    # Quadratic pair: two disjoint sign-changing rational intervals exhaust degree two.
    for lo,hi in [(F(1,5),F(1,4)),(F(3,4),F(4,5))]:check(value(q,lo)*value(q,hi)<0,'two complete rational root brackets')
    gm=[F(1),F(1,2)]
    for k in range(7):gm.append(gm[-1]-gm[-2]/6)
    check(gm[:5]==fourth_prefix+[F(7,36)],'equal-weight quadratic-root recurrence moments')
    check(interval_moments(gm[:5])['status']=='UNIQUE','lower endpoint zero-square uniqueness')
    for k in range(7):check(L([0]*k+q,gm)==0,'quadratic annihilator all supplied moments')
    # p_gamma: direct coefficient cancellation, exact negative root values, and optimum.
    for gamma in [F(-2),F(-1,17),F(0),F(1,3),F(1),F(4)]:
        p=add(lower,scale(upper,gamma))
        G0=outer(q);G1=[[gamma*x for x in row] for row in outer(ell)]
        check(gram_certificate(p,0,1,2,'even',G0,G1)==(gamma>=0),'parameter family exact Gram acceptance')
        check(L(p,gm)==gamma/F(72),'negative parameter root-pair witness')
    p1=add(lower,upper);check(p1==[F(1,36),F(-1,12),F(1,12)],'actual quadratic degree after cancellation')
    check(add(p1,[-F(1,144)])==scale(mul(ell,ell),F(1,12)) and value(p1,F(1,2))==F(1,144),'attained minimum primal point and square lower bound')
    for m in [[1,F(1,2),F(1,8)],[1,F(1,2),F(5,2)]]:
        cert=interval_moments(m);check(cert['status']=='INFEASIBLE' and L(cert['witness'],m)<0,'two terminal infeasible tables')
    bad=[F(1),F(1,2),F(1,8)]
    for k in range(3):
        for j in range(3-k):
            h=[F(1)]
            for _ in range(j):h=mul(h,[1,-1])
            check(L([0]*k+h,bad)>=0,'all available elementary differences pass false table')
    # Affine terminal migration X=2+3Y; binomial coefficients computed exactly.
    from math import comb
    low_x=[sum((F(comb(k,j))*2**(k-j)*3**j*gm[j] for j in range(k+1)),F(0)) for k in range(5)]
    high_x=atoms([2,F(7,2),5],[F(1,6),F(2,3),F(1,6)],4)
    check(low_x==[1,F(7,2),13,F(203,4),F(823,4)],'affine lower moments')
    check(high_x==[1,F(7,2),13,F(203,4),F(1655,8)],'affine upper moments')
    qx=[F(23,2),-7,1];gx=[-10,7,-1];lx=[F(-7,2),1]
    for m in [low_x,high_x]:
        check(L(mul(qx,qx),m)==m[4]-F(823,4),'affine lower dual polynomial')
        check(L(mul(gx,mul(lx,lx)),m)==F(1655,8)-m[4],'affine upper dual polynomial')
        check(interval_moments(m,2,5)['status']=='UNIQUE','affine endpoint uniqueness')
    uniform=[F(1,k+1) for k in range(5)]
    combo=[F(3,5)*gm[k]+F(2,5)*upper_m[k] for k in range(5)]
    check(combo==uniform,'same complete uniform prefix finite representative')
    # Degree-zero, zero-mass and collapsed support interfaces.
    for d in range(9):
        check(interval_moments([0]*(d+1))['status']=='ZERO','zero measure every order')
        m=atoms([F(2)],[F(3)],d);check(interval_moments(m,2,2)['status']=='UNIQUE','collapsed interval unique even at degree zero')
        if d:
            m[-1]+=1;cert=interval_moments(m,2,2);check(cert['status']=='INFEASIBLE' and L(cert['witness'],m)<0,'collapsed interval mismatch witness')
    check(interval_moments([1])['status']=='MULTIPLE','degree-zero positive mass nondegenerate interval')
    check(interval_moments([-1])['status']=='INFEASIBLE','negative mass')
    check(interval_moments([0,0,1])['status']=='INFEASIBLE','zero mass nonzero high moment')
    check(gram_certificate([0],0,1,0,'even',[[0]],[]) and gram_certificate([3],0,1,0,'even',[[3]],[]),'zero and constant Gram contracts')
    check(not rational_atom_certificate([1],0,1,[F(-1)],[1]),'out-of-support supplied atom refused')
    for fn in [lambda:interval_moments([],0,1),lambda:interval_moments([1],1,0),lambda:gram_certificate([1],0,0,0,'even',[[1]],[]),lambda:gram_certificate([0,0,1],0,1,0,'even',[[1]],[]),lambda:psd([[1,2],[0,1]])]:
        try:fn()
        except ValueError:check(True,'invalid contract refused')
        else:check(False,'invalid contract refused')
    return {'status':'PASS','checks':sum(COUNTS.values()),'groups':dict(COUNTS),'third_moment_interval':['2/9','5/18'],'fourth_moment_interval':['7/36','5/24'],'quadratic_root_brackets':[['1/5','1/4'],['3/4','4/5']],'lower_quadratic_pair_moments':[str(x) for x in gm[:5]],'upper_three_atom_moments':[str(x) for x in upper_m],'parameter_condition':'gamma >= 0','p1_minimum':'1/144','p1_minimizer':'1/2','scope':'Exact supplied Gram/moment/atom certificates and complete low-degree endpoints; no general SDP optimization or algebraic root solver.'}
if __name__=='__main__':
    ap=argparse.ArgumentParser(description=__doc__);ap.add_argument('--output',type=Path,required=True);args=ap.parse_args();result=run();args.output.parent.mkdir(parents=True,exist_ok=True);args.output.write_text(json.dumps(result,ensure_ascii=False,indent=2)+'\n');print(json.dumps(result,ensure_ascii=False))