#!/usr/bin/env python3
"""Exact Gaussian-rational certificates for finite disk and Taylor interpolation.
Standard library only; no floating eigensolver, finite mesh norm certificate, or
unproved root guesses. Repeated inputs must be merged before sample solving.
Use --output to keep independent review results out of the author directory.
"""
from dataclasses import dataclass
from fractions import Fraction as F
from itertools import product, permutations
from collections import Counter
from pathlib import Path
import argparse,json
CHECKS=Counter()
def require(c,msg):
    if not c: raise ValueError(msg)
def check(c,label):
    CHECKS[label]+=1
    if not c: raise RuntimeError(label)
def rational(x):
    require(type(x) is int or isinstance(x,F),'exact rational input required')
    return F(x)
@dataclass(frozen=True)
class C:
    re:F=F(0)
    im:F=F(0)
    def __post_init__(self):
        object.__setattr__(self,'re',rational(self.re));object.__setattr__(self,'im',rational(self.im))
    def __add__(self,b):
        b=Cv(b);return C(self.re+b.re,self.im+b.im)
    __radd__=__add__
    def __neg__(self):return C(-self.re,-self.im)
    def __sub__(self,b):return self+-Cv(b)
    def __rsub__(self,b):return Cv(b)+-self
    def __mul__(self,b):
        b=Cv(b);return C(self.re*b.re-self.im*b.im,self.re*b.im+self.im*b.re)
    __rmul__=__mul__
    def conjugate(self):return C(self.re,-self.im)
    def norm2(self):return self.re*self.re+self.im*self.im
    def __truediv__(self,b):
        b=Cv(b);q=b.norm2();require(q!=0,'zero denominator');v=self*b.conjugate();return C(v.re/q,v.im/q)
    def __rtruediv__(self,b):return Cv(b)/self
    def __bool__(self):return bool(self.re or self.im)
    def __eq__(self,b):
        try:b=Cv(b)
        except ValueError:return False
        return self.re==b.re and self.im==b.im
    def __hash__(self):return hash((self.re,self.im))
def Cv(x):return x if isinstance(x,C) else C(x)
Z=C();O=C(1);I=C(0,1)
def phi(a,z):return (z-a)/(1-a.conjugate()*z)
def trim(p):
    p=list(map(Cv,p));require(bool(p),'empty polynomial')
    while len(p)>1 and not p[-1]:p.pop()
    return p
def padd(p,q):return trim([(p[j] if j<len(p) else Z)+(q[j] if j<len(q) else Z) for j in range(max(len(p),len(q)))])
def pscale(p,c):return trim([Cv(c)*x for x in p])
def pmul(p,q):
    r=[Z]*(len(p)+len(q)-1)
    for j,a in enumerate(p):
        for k,b in enumerate(q):r[j+k]+=a*b
    return trim(r)
def pvalue(p,z):
    v=Z
    for c in reversed(p):v=v*z+c
    return v
def pdiv(p,q):
    p=trim(p);q=trim(q);require(any(q),'zero divisor');r=[Z]*max(1,len(p)-len(q)+1)
    while any(p) and len(p)>=len(q):
        j=len(p)-len(q);c=p[-1]/q[-1];r[j]+=c;p=padd(p,pscale([Z]*j+q,-c))
    return trim(r),p
def pgcd(p,q):
    while any(q):p,q=q,pdiv(p,q)[1]
    return pscale(p,1/p[-1])
def rf(n,d):
    n=trim(n);d=trim(d);require(any(d),'zero rational denominator')
    g=pgcd(n,d);n=pdiv(n,g)[0];d=pdiv(d,g)[0]
    require(bool(d[0]),'rational function must be regular at zero')
    return (pscale(n,1/d[0]),pscale(d,1/d[0]))
def rconst(c):return rf([Cv(c)],[O])
def rvalue(r,z):return pvalue(r[0],z)/pvalue(r[1],z)
def radd(r,s):return rf(padd(pmul(r[0],s[1]),pmul(s[0],r[1])),pmul(r[1],s[1]))
def rmul(r,s):return rf(pmul(r[0],s[0]),pmul(r[1],s[1]))
def rscale(r,c):return rf(pscale(r[0],c),r[1])
def rcompose(r,s):
    def poly_at(p):
        out=rconst(0)
        for c in reversed(p):out=radd(rmul(out,s),rconst(c))
        return out
    n,d=poly_at(r[0]),poly_at(r[1]);return rf(pmul(n[0],d[1]),pmul(n[1],d[0]))
def degree(r):return max(len(r[0]),len(r[1]))-1
def augment(a,b,g):
    a,b=Cv(a),Cv(b);require(a.norm2()<1 and b.norm2()<1,'strict disk augmentation parameters')
    x=pmul([-a,O],g[0]);y=pmul([O,-a.conjugate()],g[1])
    return rf(padd(pscale(y,b),x),padd(y,pscale(x,b.conjugate())))
def series(r,n):
    require(type(n) is int and n>=0,'nonnegative order');a,b=r;out=[]
    for k in range(n+1):out.append(((a[k] if k<len(a) else Z)-sum((b[j]*out[k-j] for j in range(1,min(k,len(b)-1)+1)),Z))/b[0])
    return out
def eye(n):return [[O if i==j else Z for j in range(n)] for i in range(n)]
def star(A):return [[x.conjugate() for x in row] for row in zip(*A)]
def mm(A,B):return [[sum((a*b for a,b in zip(row,col)),Z) for col in zip(*B)] for row in A]
def msub(A,B):return [[a-b for a,b in zip(x,y)] for x,y in zip(A,B)]
def mscale(A,c):return [[c*x for x in row] for row in A]
def quad(A,v):return sum((v[i].conjugate()*A[i][j]*v[j] for i in range(len(v)) for j in range(len(v))),Z)
def psd(A):
    n=len(A);require(all(len(row)==n for row in A) and A==star(A),'Hermitian square matrix required')
    if not n:return dict(psd=True,rank=0,pivots=[])
    for i in range(n):
        if A[i][i].re<0:
            v=[Z]*n;v[i]=O;return dict(psd=False,witness=v,value=A[i][i])
        if not A[i][i]:
            for j in range(n):
                if A[i][j]:
                    v=[Z]*n;v[j]=O;v[i]=-A[i][j]*((A[j][j].re+1)/(2*A[i][j].norm2()))
                    return dict(psd=False,witness=v,value=quad(A,v))
    piv=next((i for i in range(n) if A[i][i]),None)
    if piv is None:return dict(psd=True,rank=0,pivots=[])
    rem=[j for j in range(n) if j!=piv];d=A[piv][piv]
    S=[[A[i][j]-A[i][piv]*A[piv][j]/d for j in rem] for i in rem];r=psd(S)
    if r['psd']:return dict(psd=True,rank=1+r['rank'],pivots=[d]+r['pivots'])
    v=[Z]*n
    for j,x in zip(rem,r['witness']):v[j]=x
    v[piv]=-sum((A[piv][j]*v[j] for j in rem),Z)/d
    return dict(psd=False,witness=v,value=quad(A,v))
def sample_data(z,w):
    require(len(z)==len(w) and len(z)>0,'matching nonempty sample lists');z=list(map(Cv,z));w=list(map(Cv,w))
    require(all(x.norm2()<1 for x in z),'strictly interior sample nodes');require(len(set(z))==len(z),'merge repeated input nodes first')
    return z,w
def pick_matrix(z,w):
    z,w=sample_data(z,w);return [[(1-a*b.conjugate())/(1-x*y.conjugate()) for y,b in zip(z,w)] for x,a in zip(z,w)]
def sample_solve(z,w):
    z,w=sample_data(z,w);steps=[];layers=[]
    while z:
        layers.append(dict(nodes=z[:],values=w[:]))
        if any(v.norm2()>1 for v in w):return dict(status='INFEASIBLE',steps=steps,layers=layers)
        a,b=z[0],w[0]
        if b.norm2()==1:
            return dict(status='UNIQUE' if all(v==b for v in w) else 'INFEASIBLE',steps=steps,layers=layers,terminal=b)
        steps.append((a,b));w=[phi(b,v)/phi(a,x) for x,v in zip(z[1:],w[1:])];z=z[1:]
    return dict(status='FREE',steps=steps,layers=layers)
def sample_construct(result,tail=Z):
    require(result['status']!='INFEASIBLE','infeasible sample data');tail=Cv(tail);require(tail.norm2()<=1,'Schur constant tail')
    r=rconst(result['terminal'] if result['status']=='UNIQUE' else tail)
    for a,b in reversed(result['steps']):r=augment(a,b,r)
    return r
def coefficient_defect(c):
    c=list(map(Cv,c));require(bool(c),'nonempty coefficient list');n=len(c)
    T=[[c[i-j] if i>=j else Z for j in range(n)] for i in range(n)]
    return msub(eye(n),mm(T,star(T)))
def strip(c):
    g=c[0];require(g.norm2()<1,'strict coefficient pivot');d=[]
    for k in range(len(c)-1):
        d.append((c[k+1]+g.conjugate()*sum((c[j]*d[k-j] for j in range(1,k+1)),Z))/(1-g.norm2()))
    return d
def coefficient_solve(c):
    c=list(map(Cv,c));require(bool(c),'nonempty coefficient list');params=[];layers=[]
    while c:
        layers.append(c[:]);g=c[0]
        if g.norm2()>1:return dict(status='INFEASIBLE',params=params,layers=layers)
        if g.norm2()==1:return dict(status='UNIQUE' if not any(c[1:]) else 'INFEASIBLE',params=params,layers=layers,terminal=g)
        params.append(g);c=strip(c)
    return dict(status='FREE',params=params,layers=layers)
def coefficient_construct(result,tail=Z):
    require(result['status']!='INFEASIBLE','infeasible coefficient data');tail=Cv(tail);require(tail.norm2()<=1,'Schur constant tail')
    r=rconst(result['terminal'] if result['status']=='UNIQUE' else tail)
    for g in reversed(result['params']):r=augment(Z,g,r)
    return r
def blaschke(zeros,phase=O):
    zeros=list(map(Cv,zeros));phase=Cv(phase);require(phase.norm2()==1 and all(a.norm2()<1 for a in zeros),'valid Blaschke parameters')
    r=rconst(phase)
    for a in zeros:r=rmul(r,rf([-a,O],[O,-a.conjugate()]))
    return r
def derivative(r,z):
    a,b=r;ap=[k*a[k] for k in range(1,len(a))] or [Z];bp=[k*b[k] for k in range(1,len(b))] or [Z]
    return (pvalue(ap,z)*pvalue(b,z)-pvalue(a,z)*pvalue(bp,z))/(pvalue(b,z)*pvalue(b,z))
def rejected(f):
    try:f()
    except ValueError:return True
    return False

def check_sample(z,w):
    A=pick_matrix(z,w);p=psd(A);s=sample_solve(z,w)
    check(p['psd']==(s['status']!='INFEASIBLE'),'sample matrix/algorithm agreement')
    if not p['psd']:
        q=quad(A,p['witness']);check(q==p['value'] and q.im==0 and q.re<0,'negative Pick witness');return s
    check((p['rank']<len(z))==(s['status']=='UNIQUE'),'Pick uniqueness/rank')
    for tail in [Z,O,C(F(-1,2)),I]:
        f=sample_construct(s,tail);check(all(rvalue(f,Cv(x))==Cv(v) for x,v in zip(z,w)),'reconstructed all sample values')
        if s['status']=='UNIQUE':check(degree(f)==p['rank'],'unique Pick Blaschke degree')
    for layer in s['layers']:
        z1,w1=layer['nodes'],layer['values'];a,b=z1[0],w1[0]
        if b.norm2()>=1 or any(v.norm2()>1 for v in w1):continue
        if len(z1)==1:continue
        u=[phi(b,v) for v in w1];Q=pick_matrix(z1,u);t=[phi(a,x) for x in z1[1:]];eta=[v/x for v,x in zip(u[1:],t)]
        R=pick_matrix(z1[1:],eta);S=[[Q[i][j]-Q[i][0]*Q[0][j]/Q[0][0] for j in range(1,len(z1))] for i in range(1,len(z1))]
        check(S==[[t[i]*R[i][j]*t[j].conjugate() for j in range(len(t))] for i in range(len(t))],'Pick Schur complement identity')
        P=pick_matrix(z1,w1);D=[1/(1-b.conjugate()*v) for v in w1]
        check(Q==[[(1-b.norm2())*D[i]*P[i][j]*D[j].conjugate() for j in range(len(D))] for i in range(len(D))],'Pick output congruence')
    return s

def check_coeff(c):
    c=list(map(Cv,c));A=coefficient_defect(c);p=psd(A);s=coefficient_solve(c)
    check(p['psd']==(s['status']!='INFEASIBLE'),'coefficient matrix/algorithm agreement')
    if not p['psd']:
        q=quad(A,p['witness']);check(q==p['value'] and q.im==0 and q.re<0,'negative coefficient witness');return s
    check((p['rank']<len(c))==(s['status']=='UNIQUE'),'coefficient uniqueness/rank')
    for tail in [Z,O,I]:
        f=coefficient_construct(s,tail);check(series(f,len(c)-1)==c,'reconstructed all coefficients')
        if s['status']=='UNIQUE':check(degree(f)==p['rank'],'unique coefficient Blaschke degree')
    for layer in s['layers']:
        if layer[0].norm2()>=1:continue
        n=len(layer);d=strip(layer);g=layer[0]
        W=[[d[i-j-1] if i>j else Z for j in range(n)] for i in range(n)]
        E=[[eye(n)[i][j]+g.conjugate()*W[i][j] for j in range(n)] for i in range(n)]
        small=coefficient_defect(d) if d else [];block=[[Z]*n for _ in range(n)];block[0][0]=O
        for i in range(1,n):
            for j in range(1,n):block[i][j]=small[i-1][j-1]
        check(mm(mm(E,coefficient_defect(layer)),star(E))==mscale(block,1-g.norm2()),'finite Toeplitz congruence')
    return s

def run():
    CHECKS.clear();nodes=list(map(Cv,[F(-1,2),0,F(1,2)]));records={}
    for t in [F(1,3),F(2,5),F(9,20)]:
        w=[Z,Z,C(t)];s=check_sample(nodes,w);records[str(t)]={'matrix':pick_matrix(nodes,w),'decision':s['status'],'psd':psd(pick_matrix(nodes,w))}
    P=pick_matrix(nodes,[Z,Z,C(F(9,20))]);v=list(map(Cv,[F(3,5),F(-8,5),1]))
    check(quad(P,v)==C(F(-17,300)),'terminal three-point obstruction')
    for a,b in [(0,1),(0,2),(1,2)]:check(psd([[P[i][j] for j in [a,b]] for i in [a,b]])['psd'],'all pairs still feasible')
    B=blaschke([Z,C(F(-1,2))]);check(rvalue(B,C(F(1,2)))==C(F(2,5)),'two-zero exact radius')
    f0=rscale(B,F(5,6));f1=rmul(B,rf([4,7],[7,4]))
    check(rvalue(f0,C(F(1,4)))==C(F(5,36)) and rvalue(f1,C(F(1,4)))==C(F(23,192)),'distinct positive Pick solutions')
    check(rvalue(rscale(B,F(9,8)),C(F(1,2)))==C(F(9,20)),'minimum norm attaining function')
    # Same conclusion after complex input/output automorphisms.
    a,b=C(F(1,4),F(1,3)),C(F(-1,3),F(1,5));zn=[phi(a,z) for z in nodes]
    for t in [F(1,3),F(2,5),F(9,20)]:
        wn=[phi(b,v) for v in [Z,Z,C(t)]];s=check_sample(zn,wn);check(s['status']==records[str(t)]['decision'],'complex coordinate migration')
    for delta,status in [(F(-1,1000),'FREE'),(F(1,1000),'INFEASIBLE')]:
        t=F(2,5)+delta;s=check_sample(nodes,[Z,Z,C(t)]);check(s['status']==status,'near-boundary exact classification')
    psi=rf([-I,O],[I,O]);upper_nodes=[I/3,I,3*I]
    check([rvalue(psi,x) for x in upper_nodes]==nodes,'half-plane node migration')
    bh=rcompose(B,psi);explicit=rf(pmul([-I,O],[-I,C(3)]),pmul([I,O],[I,C(3)]))
    check(bh==explicit and derivative(rscale(bh,F(5,6)),I)==C(0,F(-5,24)),'half-plane rational and derivative')
    coef_records={}
    for c2 in [F(3,32),F(15,32),F(1,2)]:
        c=list(map(Cv,[F(1,2),F(3,8),c2]));s=check_coeff(c);coef_records[str(c2)]={'decision':s['status'],'psd':psd(coefficient_defect(c))}
    strict=coefficient_construct(coefficient_solve([F(1,2),F(3,8),F(3,32)]));unique=coefficient_construct(coefficient_solve([F(1,2),F(3,8),F(15,32)]))
    check(strict==rf([6,7,4],[12,5,2]),'strict rational endpoint')
    check(unique==rf([2,3,4],[4,3,2]) and degree(unique)==2,'unique rational endpoint')
    check(pvalue(list(map(Cv,[F(1,2),F(3,8),F(15,32)])),C(F(4,5)))==C(F(11,10)) and rvalue(unique,C(F(4,5)))==C(F(29,32)),'truncation failure and true extension')
    # Exhaustive finite coefficient box includes nonreal parameters and degeneracies.
    grid=list(map(Cv,[F(-1),F(-1,2),0,F(1,2),1]))+[I,-I,C(F(1,3),F(1,4))]
    for c in product(grid,repeat=3):check_coeff(c)
    # Exhaustive real three-output sample box, with an independent Hermitian test.
    for w in product([C(F(k,4)) for k in range(-4,5)],repeat=3):check_sample(nodes,w)
    family_nodes=[Z,C(F(1,2)),C(F(-1,3)),C(0,F(1,2))]
    families=[rconst(C(F(1,3),F(1,4))),blaschke([Z]),blaschke([C(F(1,3),F(1,4)),Z]),rscale(blaschke([Z,C(F(-1,2))]),F(2,3)),unique,strict]
    for f in families:
        for perm in permutations(family_nodes):check_sample(list(perm),[rvalue(f,z) for z in perm])
        for order in range(6):check_coeff(series(f,order))
    circle=[C((1-F(k,3)**2)/(1+F(k,3)**2),2*F(k,3)/(1+F(k,3)**2)) for k in range(-12,13)]+[-O]
    for zeros in [[Z],[Z,Z],[Z,C(F(1,2))],[C(F(1,3),F(1,4)),C(F(-1,3),F(1,4)),Z]]:
        f=blaschke(zeros,I);check(degree(f)==len(zeros),'Blaschke exact degree')
        for z in circle:
            value=rvalue(f,z);check(value.norm2()==1,'exact circle modulus')
            speed=z*derivative(f,z)/value;kernel=sum(((1-a.norm2())/(z-a).norm2() for a in zeros),F())
            check(speed==C(kernel) and kernel>0,'positive boundary phase derivative')
        for a in zeros:check(rvalue(f,a)==0,'specified zeros')
    two=blaschke([Z,C(F(1,2))]);check(rvalue(two,I)==C(F(-3,5),F(-4,5)),'twofold explicit circle image')
    check(derivative(two,O)==4 and (-O)*derivative(two,-O)/rvalue(two,-O)==F(4,3),'twofold phase endpoint speeds')
    check(degree(rmul(two,blaschke([Z,Z,Z])))==5 and degree(rcompose(two,blaschke([Z,Z,Z])))==6,'product versus composition degree')
    critical=list(map(Cv,[1,-4,1]));num,den=two
    num_prime=[k*num[k] for k in range(1,len(num))];den_prime=[k*den[k] for k in range(1,len(den))]
    derivative_numerator=padd(pmul(num_prime,den),pscale(pmul(num,den_prime),-1))
    check(not any(pdiv(derivative_numerator,critical)[1]),'critical derivative vanishes')
    value_numerator=padd(num,pmul(list(map(Cv,[0,0,1])),den))
    check(not any(pdiv(value_numerator,critical)[1]),'actual critical value B(c)+c squared')
    # For w=-c^2, four times the quadratic discriminant is (c^2+1)^2-16c^2.
    discriminant_numerator=padd(pmul(list(map(Cv,[1,0,1])),list(map(Cv,[1,0,1]))),list(map(Cv,[0,0,-16])))
    check(discriminant_numerator==list(map(Cv,[1,0,-14,0,1])) and not any(pdiv(discriminant_numerator,critical)[1]),'actual double preimage discriminant')
    for c2 in [F(3,32),F(15,32),F(1,2)]:
        old=coefficient_solve([F(1,2),F(3,8),c2]);rotated=[I*C(F(1,2)),-C(F(3,8)),-I*c2]
        new=check_coeff(rotated);check(new['status']==old['status'],'complex coefficient rotation')
    check_coeff([1,0,0]);check_coeff([1,0,F(1,1000)]);check_sample([0,F(1,2)],[I,I]);check_sample([0,F(1,2)],[I,0])
    for fn in [lambda:C(0.5),lambda:sample_solve([1],[0]),lambda:sample_solve([0,0],[0,0]),lambda:sample_solve([0,0],[0,1]),lambda:sample_solve([],[]),lambda:coefficient_solve([]),lambda:blaschke([1]),lambda:blaschke([0],F(1,2)),lambda:sample_construct(sample_solve([0],[2])),lambda:coefficient_construct(coefficient_solve([2]))]:
        check(rejected(fn),'invalid or infeasible input rejection')
    return dict(status='PASS',checks=sum(CHECKS.values()),groups=dict(CHECKS),sample_endpoints=records,coefficient_endpoints=coef_records,strict_function=strict,unique_function=unique,negative_pick_vector=v,scope='Exact finite checks and explicit reconstructions; the general analytic proof is in the articles, not inferred from sampled values')
def serial(x):
    if isinstance(x,C):return {'re':str(x.re),'im':str(x.im)}
    if isinstance(x,F):return str(x)
    if isinstance(x,(list,tuple)):return [serial(v) for v in x]
    if isinstance(x,dict):return {k:serial(v) for k,v in x.items()}
    return x
if __name__=='__main__':
    ap=argparse.ArgumentParser(description=__doc__);ap.add_argument('--output',type=Path);args=ap.parse_args();data=json.dumps(serial(run()),ensure_ascii=False,indent=2)+'\n'
    if args.output:args.output.parent.mkdir(parents=True,exist_ok=True);args.output.write_text(data);print('PASS',sum(CHECKS.values()),'checks;',args.output)
    else:print(data,end='')