"""Exact sextic identification and construction verification (stdlib only).
Run normally for algebra + all matching flows; --algebra-only skips the latter.
Polynomial coefficients are in ascending order. No floating point is used.
"""
from fractions import Fraction as F
from pathlib import Path
from math import gcd
from functools import reduce
import json,runpy,sys
HERE=Path(__file__).parent

def trim(p):
 p=list(p)
 while len(p)>1 and p[-1]==0:p.pop()
 return p or [F(0)]
def add(a,b):return trim([(a[i]if i<len(a)else 0)+(b[i]if i<len(b)else 0)for i in range(max(len(a),len(b)))])
def neg(a):return [-v for v in a]
def sub(a,b):return add(a,neg(b))
def mul(a,b):
 r=[F(0)]*(len(a)+len(b)-1)
 for i,x in enumerate(a):
  for j,y in enumerate(b):r[i+j]+=x*y
 return trim(r)
def divmodp(a,b):
 a=trim(list(map(F,a)));b=trim(list(map(F,b)));assert b!=[0]
 q=[F(0)]*max(1,len(a)-len(b)+1)
 while a!=[0]and len(a)>=len(b):
  k=len(a)-len(b);v=a[-1]/b[-1];q[k]=v
  a=sub(a,[F(0)]*k+[v*x for x in b])
 return trim(q),a
def derivative(a):return trim([i*a[i]for i in range(1,len(a))])
def evaluate(p,x):
 v=0
 for c in reversed(p):v=v*x+c
 return v
def sturm_count(p,lo,hi):
 seq=[list(map(F,p)),derivative(list(map(F,p)))]
 while seq[-1]!=[0]:
  r=neg(divmodp(seq[-2],seq[-1])[1])
  if r==[0]:break
  seq.append(r)
 def variations(x):
  v=[evaluate(p,x)for p in seq];v=[1 if y>0 else -1 for y in v if y]
  return sum(a!=b for a,b in zip(v,v[1:]))
 assert evaluate(p,lo)and evaluate(p,hi)
 return variations(lo)-variations(hi)

PT=list(map(F,[105476,-423516,422009,-226450,211627,-148176,34006]))
PG=list(map(F,[20566231563858477973,-315750351353916573206,2019436570431176007853,-6886878440990980289832,13208173323262094551248,-13507216477071755171328,5754135102948532776960]))
DP=list(map(F,[285,469,-1512,694]))
KP=list(map(F,[16,77,-220,111]))
GN=list(map(F,[1484,2072,-7130,3367,-49]))
AN=mul([-1,1],[-14,-32,-6,7])
NUMS=[
 [226,-389,312,-77],
 [637,675,-2981,1469],
 [452,1541,-3330,1337],
 [32,154,-440,222],
 [250,886,-1995,835],
 [113,-53,-638,378]]

# Rational function field, with t subject to PT(t)=0.
class E:
 def __init__(self,x):
  if isinstance(x,E):self.p=x.p;return
  self.p=divmodp(x if isinstance(x,list)else [F(x)],PT)[1]
 def __add__(self,x):return E(add(self.p,E(x).p))
 __radd__=__add__
 def __neg__(self):return E(neg(self.p))
 def __sub__(self,x):return self+-E(x)
 def __rsub__(self,x):return E(x)+-self
 def __mul__(self,x):return E(mul(self.p,E(x).p))
 __rmul__=__mul__
 def inverse(self):
  oldr,r=PT,self.p;olds,ss=[F(0)],[F(1)]
  while r!=[0]:
   q,rem=divmodp(oldr,r);oldr,r=r,rem;olds,ss=ss,sub(olds,mul(q,ss))
  assert len(oldr)==1 and oldr[0]
  return E([x/oldr[0]for x in olds])
 def __truediv__(self,x):return self*E(x).inverse()
 def __rtruediv__(self,x):return E(x)*self.inverse()
 def __eq__(self,x):return self.p==E(x).p

t=E([0,1]);D=E(DP);K=E(KP);alpha=E(GN)/(12*D)
a=[E(p)/(6*D)for p in NUMS]
lam=[E(AN)/(2*K),t];lam.append(1-lam[0]-lam[1]);T=1-alpha
assert sum(a)==1 and sum(lam)==1
assert sub(mul(derivative(GN),DP),mul(GN,derivative(DP)))==neg(PT)
assert evaluate(PG,alpha)==0
cert=json.loads((HERE/'stationarity-input.json').read_text())
Cs=[[[F(x,4)for x in row]for row in C]for C in cert['C4']]
Qs=[sum(a[i]*C[i][j]*a[j]for i in range(6)for j in range(6))for C in Cs]
assert all(Q==T for Q in Qs)
for i in range(6):assert sum(lam[k]*Cs[k][i][j]*a[j]for k in range(3)for j in range(6))==T
print('PASS: exact mass normalization, three active values, and stationarity modulo P(t).')
print('PASS: exact derivative identity and minimal-polynomial composition.')

# Irreducibility over F_41 by Rabin's criterion (degree n=6).
def mpdiv(a,b,p):
 a=trim([int(x)%p for x in a]);b=trim([int(x)%p for x in b]);q=[0]*max(1,len(a)-len(b)+1)
 while a!=[0]and len(a)>=len(b):
  k=len(a)-len(b);v=a[-1]*pow(b[-1],-1,p)%p;q[k]=v
  for j in range(len(b)):a[j+k]=(a[j+k]-v*b[j])%p
  a=trim(a)
 return trim(q),a
def mppow(a,k,mod,p):
 r=[1]
 while k:
  if k&1:r=mpdiv(mul(r,a),mod,p)[1]
  a=mpdiv(mul(a,a),mod,p)[1];k//=2
 return r
def mpgcd(a,b,p):
 while b!=[0]:a,b=b,mpdiv(a,b,p)[1]
 return a
for P in(PT,PG):
 assert reduce(gcd,[abs(int(x))for x in P])==1 and int(P[-1])%41
 pp=mpdiv(P,[1],41)[0] # coefficient reduction, retaining the full polynomial
 assert len(pp)==7
 assert mpdiv(sub(mppow([0,1],41**6,pp,41),[0,1]),pp,41)[1]==[0]
 for k in(2,3):assert len(mpgcd(pp,mpdiv(sub(mppow([0,1],41**k,pp,41),[0,1]),pp,41)[1],41))==1
print('PASS: both primitive sextics are irreducible modulo 41, hence over Q.')

class I:
 def __init__(self,l,h=None):self.l=F(l);self.h=F(l if h is None else h)
 def __add__(self,x):x=ii(x);return I(self.l+x.l,self.h+x.h)
 __radd__=__add__
 def __neg__(self):return I(-self.h,-self.l)
 def __sub__(self,x):return self+-ii(x)
 def __rsub__(self,x):return ii(x)+-self
 def __mul__(self,x):
  x=ii(x);v=[self.l*x.l,self.l*x.h,self.h*x.l,self.h*x.h];return I(min(v),max(v))
 __rmul__=__mul__
 def __truediv__(self,x):
  x=ii(x);assert x.l*x.h>0;return self*I(1/x.h,1/x.l)
 def __rtruediv__(self,x):return ii(x)/self

def ii(x):return x if isinstance(x,I)else I(x)
# Decimal strings here are exact rational isolating bounds, not floating inputs.
root_data=json.loads((HERE/'root-intervals.json').read_text());tl,th=map(F,root_data['t_interval']);gl,gh=map(F,root_data['alpha_interval'])
assert sturm_count(PT,tl,th)==1
assert sturm_count(PT,F('0.359275945611134729'),F('0.359275945611134730'))==1
assert sturm_count(PG,gl,gh)==1
assert sturm_count(PG,F('0.41977282048624007772'),F('0.41977282048624007773'))==1
TI=I(tl,th);DI=evaluate(DP,TI);KI=evaluate(KP,TI)
assert DI.l>0 and KI.l>0
ai=[evaluate(p,TI)/(6*DI)for p in NUMS]
li=[evaluate(AN,TI)/(2*KI),TI];li.append(1-li[0]-li[1])
gI=evaluate(GN,TI)/(12*DI)
assert gl<gI.l<=gI.h<gh
assert all(x.l>0 for x in ai+li)
point=ai+li[:2]+[1-gI];center=list(map(F,cert['center']));radius=F(cert['radius'])
for val,c in zip(point,center):assert c-radius<val.l<=val.h<c+radius
print('PASS: Sturm root isolation, positivity, and inclusion in the certified flow box.')
print('Exact alpha lies between',root_data['alpha_interval'][0],'and',root_data['alpha_interval'][1])
if '--algebra-only' not in sys.argv:
 runpy.run_path(str(HERE/'verify_flow_bound.py'),run_name='__main__')
 print('PASS: gamma(4) >= alpha for the exact algebraic six-block construction.')