#!/usr/bin/env python3
"""Exact Gaussian-rational contour and coefficient-family certificates.
The implementation uses no floating point, external solver, atan, or root guesses.
It implements a sufficient Pellet test and elementary tube subdivision, not the
accelerated CIsolate algorithm cited in the accompanying theory.
"""
from dataclasses import dataclass
from fractions import Fraction as F
from math import comb
from pathlib import Path
from collections import Counter
from itertools import product
import argparse,copy,json
CHECKS=Counter()
def require(v,message):
 if not v:raise ValueError(message)
def check(v,label):
 CHECKS[label]+=1
 if not v:raise RuntimeError(label)
def rational(v):
 require(type(v) is int or isinstance(v,F),'exact rational required')
 return F(v)
@dataclass(frozen=True)
class G:
 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,v):
  v=gauss(v);return G(self.re+v.re,self.im+v.im)
 __radd__=__add__
 def __neg__(self):return G(-self.re,-self.im)
 def __sub__(self,v):return self+-gauss(v)
 def __rsub__(self,v):return gauss(v)+-self
 def __mul__(self,v):
  v=gauss(v);return G(self.re*v.re-self.im*v.im,self.re*v.im+self.im*v.re)
 __rmul__=__mul__
 def __pow__(self,n):
  require(type(n) is int and n>=0,'nonnegative integer power required');a=self;b=G(1)
  while n:
   if n%2:b=b*a
   a=a*a;n//=2
  return b
 def __bool__(self):return bool(self.re or self.im)
def gauss(v):return v if isinstance(v,G) else G(rational(v))
def upper(z):z=gauss(z);return abs(z.re)+abs(z.im)
def lower(z):z=gauss(z);return max(abs(z.re),abs(z.im))
def norm2(z):z=gauss(z);return z.re*z.re+z.im*z.im
def cross(a,b):return a.re*b.im-a.im*b.re
def dot(a,b):return a.re*b.re+a.im*b.im
def polynomial(p):
 require(isinstance(p,(list,tuple)) and len(p)>0,'nonempty coefficients required');a=list(map(gauss,p))
 while len(a)>1 and not a[-1]:a.pop()
 require(any(a),'zero polynomial excluded');return tuple(a)
def value(p,z):
 z=gauss(z);v=G()
 for a in reversed(p):v=v*z+a
 return v
def multiply(p,q):
 out=[G()]*(len(p)+len(q)-1)
 for j,a in enumerate(p):
  for k,b in enumerate(q):out[j+k]=out[j+k]+a*b
 return tuple(out)
def shift(p,c):
 # Polynomial Horner in w+c: O(n^2) Gaussian arithmetic, O(n) workspace.
 c=gauss(c);out=[p[-1]]
 for a in reversed(p[:-1]):
  nxt=[G()]*(len(out)+1)
  for j,b in enumerate(out):nxt[j]=nxt[j]+c*b;nxt[j+1]=nxt[j+1]+b
  nxt[0]=nxt[0]+a;out=nxt
 return tuple(out)
def rectangle(x0,x1,y0,y1):
 x0,x1,y0,y1=map(rational,(x0,x1,y0,y1));require(x0<x1 and y0<y1,'positive rectangle required')
 return (G(x0,y0),G(x1,y0),G(x1,y1),G(x0,y1))
def curve_check(vertices):
 require(isinstance(vertices,(list,tuple)) and len(vertices)>=1,'nonempty closed polygon required')
 return tuple(map(gauss,vertices))
def polygon_winding(vertices):
 vs=curve_check(vertices);w=0
 for a,b in zip(vs,vs[1:]+vs[:1]):
  require(bool(a) and bool(b),'image vertex is zero')
  d=cross(a,b)
  require(not(d==0 and dot(a,b)<=0),'image segment meets zero')
  if a.im<=0<b.im and d>0:w+=1
  elif b.im<=0<a.im and d<0:w-=1
 return w

def segment_data(p,a,b):
 length=upper(b-a);cs=shift(p,a)
 derivative=sum((F(j)*upper(cs[j])*length**(j-1) for j in range(1,len(cs))),F())
 radius=length*derivative;floor=lower(cs[0])
 return dict(start_value=cs[0],end_value=value(p,b),length_upper=length,derivative_upper=derivative,image_radius=radius,modulus_lower=floor,margin=floor-radius)
def contour_certificate(p,vertices,max_depth=12):
 p=polynomial(p);vs=curve_check(vertices);require(type(max_depth) is int and max_depth>=0,'nonnegative integer depth required')
 leaves=[];tested=0
 for edge in range(len(vs)):
  start=vs[edge];step=vs[(edge+1)%len(vs)]-start;stack=[(F(),F(1),0)]
  while stack:
   t0,t1,depth=stack.pop();a=start+step*t0;b=start+step*t1;data=segment_data(p,a,b);tested+=1
   if not data['start_value'] or not data['end_value']:status='BOUNDARY_ZERO'
   elif data['margin']>0:status='CERTIFIED'
   elif depth>=max_depth:status='UNRESOLVED'
   else:
    mid=(t0+t1)/2;stack.extend([(mid,t1,depth+1),(t0,mid,depth+1)]);continue
   leaves.append(dict(edge=edge,t0=t0,t1=t1,depth=depth,status=status,**data))
 result=dict(status='CERTIFIED' if all(r['status']=='CERTIFIED' for r in leaves) else 'UNRESOLVED',tested_segments=tested,leaves=leaves)
 result['winding']=polygon_winding(tuple(r['start_value'] for r in leaves)) if result['status']=='CERTIFIED' else None
 return result

def verify_contour(p,vertices,cert):
 p=polynomial(p);vs=curve_check(vertices);rows=cert['leaves'];require(bool(rows),'empty certificate')
 batches=[[] for _ in vs];previous_edge=-1
 for row in rows:
  edge=row['edge'];require(type(edge) is int and 0<=edge<len(vs),'invalid original edge')
  require(edge>=previous_edge,'invalid edge order');batches[edge].append(row);previous_edge=edge
 for edge,batch in enumerate(batches):
  require(bool(batch),'missing original edge');t=F()
  for row in batch:
   require(row['t0']==t and t<row['t1']<=1,'partition gap, overlap, or reversed interval')
   require(type(row['depth']) is int and row['depth']>=0 and row['t1']-row['t0']==F(1,2**row['depth']),'invalid dyadic depth')
   a=vs[edge]+(vs[(edge+1)%len(vs)]-vs[edge])*row['t0'];b=vs[edge]+(vs[(edge+1)%len(vs)]-vs[edge])*row['t1']
   data=segment_data(p,a,b)
   for name,x in data.items():require(row[name]==x,'altered segment '+name)
   require(row['status']=='CERTIFIED' and data['margin']>0,'uncertified leaf')
   t=row['t1']
  require(t==1,'edge not completely covered')
 require(cert['status']=='CERTIFIED','overall status not certified')
 w=polygon_winding(tuple(row['start_value'] for row in rows));require(cert['winding']==w,'incorrect winding')
 return w

def pellet_certificate(p,center,radius,k,errors=None):
 p=polynomial(p);c=gauss(center);r=rational(radius);n=len(p)-1
 require(r>0 and type(k) is int and 0<=k<=n,'invalid radius or root count')
 eps=tuple(F() for _ in p) if errors is None else tuple(map(rational,errors))
 require(len(eps)==len(p) and all(e>=0 for e in eps),'invalid coefficient errors')
 cs=shift(p,c);L=lower(cs[k])*r**k;U=sum((upper(a)*r**j for j,a in enumerate(cs) if j!=k),F())
 perturbation=sum((e*(upper(c)+r)**j for j,e in enumerate(eps)),F())
 # Upper(center) is a certified modulus bound; conservative for oblique centers.
 margin=L-U-perturbation
 return dict(center=c,radius=r,count=k,shifted_coefficients=cs,dominant_lower=L,other_upper=U,coefficient_error=perturbation,margin=margin,status='CERTIFIED' if margin>0 else 'UNRESOLVED')
def verify_clusters(p,errors,clusters):
 p=polynomial(p);eps=tuple(map(rational,errors));require(len(eps)==len(p) and all(e>=0 for e in eps),'invalid coefficient errors')
 require(lower(p[-1])>eps[-1],'degree not fixed across coefficient family')
 total=0
 for j,row in enumerate(clusters):
  rebuilt=pellet_certificate(p,row['center'],row['radius'],row['count'],eps)
  require(rebuilt==row and row['margin']>0,'incorrect or unsuccessful disk certificate')
  require(row['count']>0,'empty disks omitted from complete cluster list');total+=row['count']
  for prev in clusters[:j]:require(norm2(row['center']-prev['center'])>(row['radius']+prev['radius'])**2,'cluster closures overlap')
 require(total==len(p)-1,'root counts do not close degree')
 return total

def encode(x):
 if isinstance(x,F):return str(x)
 if isinstance(x,G):return [str(x.re),str(x.im)]
 if isinstance(x,dict):return {str(k):encode(v) for k,v in x.items()}
 if isinstance(x,(tuple,list)):return [encode(v) for v in x]
 return x
def rejected(fn,*args):
 try:fn(*args)
 except (ValueError,KeyError,TypeError):return True
 return False

def run():
 CHECKS.clear()
 # Independent binomial formula checks the Horner shift on new complex centers.
 for n in range(7):
  for phase in range(9):
   pp=tuple(G(F((j+phase)%5-2,3),F((2*j+phase)%7-3,5)) for j in range(n))+ (G(1),)
   for center in [G(),G(1),G(-1),G(1,1),G(F(-2,3),F(1,4))]:
    shifted=shift(pp,center)
    expected=tuple(sum((pp[m]*comb(m,j)*center**(m-j) for m in range(j,n+1)),G()) for j in range(n+1))
    check(shifted==expected,'horner_shift_vs_binomial_formula')
    for offset in [G(F(1,3),F(2,5)),G(-2,1)]:check(value(shifted,offset)==value(pp,center+offset),'shifted_polynomial_identity')
 diamond=(G(1),G(0,1),G(-1),G(0,-1))
 check(polygon_winding(diamond)==1 and polygon_winding(tuple(reversed(diamond)))==-1,'half_open_ray_vertices')
 check(polygon_winding((G(1),G(2),G(2,1),G(1,1)))==0,'horizontal_positive_ray_edge')
 square=rectangle(-1,1,-1,1);p4=polynomial([0,0,0,0,1])
 coarse=tuple(value(p4,z) for z in square)
 check(coarse==(G(-4),)*4 and polygon_winding(coarse)==0,'coarse_sampling_misses_four_turns')
 c4=contour_certificate(p4,square);check(verify_contour(p4,square,c4)==4,'certified_four_turns')
 reverse=contour_certificate(p4,tuple(reversed(square)));check(verify_contour(p4,tuple(reversed(square)),reverse)==-4,'orientation_reversal')
 twice=contour_certificate(p4,square+square);check(verify_contour(p4,square+square,twice)==8,'multiple_traversal')
 exhausted=contour_certificate(p4,square,0);check(exhausted['status']=='UNRESOLVED' and exhausted['winding'] is None,'budget_preserves_unknown')
 boundary=contour_certificate([-1,1],square,5)
 check(boundary['status']=='UNRESOLVED' and any(r['status']=='BOUNDARY_ZERO' for r in boundary['leaves']),'actual_boundary_zero_not_half_root')
 check(rejected(verify_contour,[-1,1],square,boundary),'boundary_certificate_rejected')
 unsampled=contour_certificate([-G(1,F(1,3)),G(1)],square,7)
 check(unsampled['status']=='UNRESOLVED' and not any(r['status']=='BOUNDARY_ZERO' for r in unsampled['leaves']),'unsampled_boundary_root_remains_unknown')
 # The root counter receives expanded coefficients; factor lists only form an oracle.
 models=0
 roots_pool=(G(),G(F(1,3)),G(0,F(-1,2)),G(2),G(-2,1))
 for length in range(1,4):
  for indices in product(range(len(roots_pool)),repeat=length):
   roots=[roots_pool[j] for j in indices];p=(G(1),)
   for root in roots:p=multiply(p,(-root,G(1)))
   cert=contour_certificate(p,square,12);expected=sum(-1<r.re<1 and -1<r.im<1 for r in roots)
   check(verify_contour(p,square,cert)==expected,'expanded_polynomial_count_vs_known_roots');models+=1
 # Main six-degree cluster family.
 p=(G(1),)
 for root in (G(1),G(1),G(1),G(-1),G(-1),G(0,2)):p=multiply(p,(-root,G(1)))
 check(p==(G(0,2),G(-1,-2),G(1,-4),G(2,4),G(-2,2),G(-1,-2),G(1)),'six_degree_expansion')
 eps=(F(1,1000),)*7
 clusters=[pellet_certificate(p,c,F(1,4),k,eps) for c,k in ((G(1),3),(G(-1),2),(G(0,2),1))]
 check(verify_clusters(p,eps,clusters)==6,'entire_coefficient_family_complete')
 check(verify_clusters([1],[F(1,2)],[])==0,'constant_family_empty_cluster_list')
 # Rational concrete perturbations cross-check counts in surrounding squares; Pellet certifies the disks.
 perturbations=[]
 for phase in range(6):
  pp=tuple(a+G(F((j+phase)%3-1,2000),F((2*j+phase)%3-1,2000)) for j,a in enumerate(p))
  check(all(norm2(a-b)<=e*e for a,b,e in zip(pp,p,eps)),'actual_perturbation_inside_family')
  counts=[]
  for c,k in ((G(1),3),(G(-1),2),(G(0,2),1)):
   box=rectangle(c.re-F(1,4),c.re+F(1,4),c.im-F(1,4),c.im+F(1,4));cc=contour_certificate(pp,box,14)
   counts.append(verify_contour(pp,box,cc));check(counts[-1]==k,'perturbed_rectangle_counts')
  perturbations.append(dict(coefficients=pp,counts=counts))
 # Meromorphic count: two zeros and two poles give net zero, not no zeros.
 numerator=multiply((-G(F(1,2)),G(1)),(-G(F(1,2)),G(1)))
 denominator=multiply((G(F(1,2)),G(1)),(G(F(1,2)),G(1)))
 cn=contour_certificate(numerator,square);cd=contour_certificate(denominator,square)
 check(verify_contour(numerator,square,cn)==2 and verify_contour(denominator,square,cd)==2,'meromorphic_zero_minus_pole')
 # Pellet failure is not absence: z^2-1 at radius 1 has boundary roots.
 fail=pellet_certificate([-1,0,1],0,1,2)
 check(fail['margin']==0 and fail['status']=='UNRESOLVED','strict_pellet_boundary')
 # Every mutation is checked from the original polynomial and contour.
 mutations=[]
 for label,mutate in [
 ('missing_leaf',lambda c:c['leaves'].pop()),
 ('altered_margin',lambda c:c['leaves'][0].__setitem__('margin',F(999))),
 ('altered_value',lambda c:c['leaves'][0].__setitem__('start_value',G(999))),
 ('altered_winding',lambda c:c.__setitem__('winding',0)),
 ('invalid_edge',lambda c:c['leaves'][0].__setitem__('edge',-1)),
 ('edge_order',lambda c:c['leaves'].reverse()),
 ('gap',lambda c:c['leaves'][0].__setitem__('t1',F())),
 ('unresolved_marked_safe',lambda c:c['leaves'][0].__setitem__('status','UNRESOLVED'))]:
  bad=copy.deepcopy(c4);mutate(bad);ok=rejected(verify_contour,p4,square,bad);check(ok,'damaged_contour_rejected');mutations.append(label)
 for label,mutate in [('missing_cluster',lambda c:c.pop()),('wrong_count',lambda c:c[0].__setitem__('count',2)),('wrong_shift',lambda c:c[0].__setitem__('shifted_coefficients',())),('wrong_error',lambda c:c[0].__setitem__('coefficient_error',F()))]:
  bad=copy.deepcopy(clusters);mutate(bad);check(rejected(verify_clusters,p,eps,bad),'damaged_cluster_rejected');mutations.append(label)
 for fn,args in [(G,(True,)),(G,(.5,)),(polynomial,([0,0],)),(rectangle,(0,0,0,1)),(polygon_winding,((G(-1),G(1)),)),(contour_certificate,(p4,square,True)),(pellet_certificate,(p,0,0,1)),(pellet_certificate,(p,0,1,True)),(verify_clusters,(p,(F(1),)*7,clusters)),(verify_clusters,([1],[F(-1)],[]))]:
  check(rejected(fn,*args),'invalid_contract_rejected')
 return encode(dict(status='PASS',checks=sum(CHECKS.values()),groups=dict(sorted(CHECKS.items())),known_root_models=models,coarse_image=coarse,main_winding=c4,reverse_winding=reverse['winding'],twice_winding=twice['winding'],budget_result=exhausted,boundary_result=boundary,cluster_polynomial=p,coefficient_errors=eps,clusters=clusters,concrete_perturbations=perturbations,meromorphic_counts=dict(zeros=2,poles=2,net=0),mutation_tests=mutations,scope='Exact finite certificate checks, not a black-box complete root isolator; no floating arithmetic, automatic precision oracle, or imported solver. General proofs are in the accompanying pages.'))
if __name__=='__main__':
 ap=argparse.ArgumentParser();ap.add_argument('--output',type=Path,required=True);args=ap.parse_args();out=run();args.output.parent.mkdir(parents=True,exist_ok=True);args.output.write_text(json.dumps(out,ensure_ascii=False,indent=2)+'\n');print(out['status'],out['checks'],'checks')