"""Independent exact audit of the elimination identities. Standard library only."""
from fractions import Fraction as Q
from itertools import permutations
from pathlib import Path
import json

class Poly:
 def __init__(self,value=0):
  self.d={m:Q(c) for m,c in value.items() if c} if isinstance(value,dict) else ({(0,0,0):Q(value)} if value else {})
 def __add__(self,other):
  other=other if isinstance(other,Poly) else Poly(other);d=self.d.copy()
  for m,c in other.d.items():d[m]=d.get(m,0)+c
  return Poly(d)
 __radd__=__add__
 def __neg__(self):return Poly({m:-c for m,c in self.d.items()})
 def __sub__(self,o):return self+-coerce(o)
 def __rsub__(self,o):return coerce(o)+-self
 def __mul__(self,o):
  o=coerce(o);d={}
  for m,c in self.d.items():
   for n,e in o.d.items():
    k=tuple(x+y for x,y in zip(m,n));d[k]=d.get(k,0)+c*e
  return Poly(d)
 __rmul__=__mul__
 def __pow__(self,n):
  z=Poly(1)
  for _ in range(n):z=z*self
  return z
 def __truediv__(self,n):return self*Q(1,n)
 def deriv(self,k):return Poly({tuple(x-int(i==k) for i,x in enumerate(m)):c*m[k] for m,c in self.d.items() if m[k]})
 def coefficient(self,k,n):return Poly({tuple(0 if i==k else x for i,x in enumerate(m)):c for m,c in self.d.items() if m[k]==n})
 def iszero(self):return not self.d
 def ev(self,values):return sum(c*product(values[i]**m[i] for i in range(3))for m,c in self.d.items())
def coerce(o):return o if isinstance(o,Poly)else Poly(o)
def product(xs):
 p=1
 for x in xs:p=p*x
 return p
def det(A):
 n=len(A);out=Poly()
 for p in permutations(range(n)):
  sign=(-1)**sum(p[i]>p[j]for i in range(n)for j in range(i+1,n))
  out+=sign*product(A[i][p[i]]for i in range(n))
 return out
u=Poly({(1,0,0):1});v=Poly({(0,1,0):1});g=Poly({(0,0,1):1})
DATA=Path(__file__).resolve().parent/'stationarity-input.json'
D=json.loads(DATA.read_text())
C=[[[Q(x,4)for x in row]for row in m]for m in D['C4']]
# Embed (p,b,z,q,f) by d=(-p+b-z+q-f)/4,e=q-d.
B=[[Q(int(i==j))for j in range(5)]for i in range(3)]
B += [[Q(-1,4),Q(1,4),Q(-1,4),Q(1,4),Q(-1,4)], [Q(1,4),Q(-1,4),Q(1,4),Q(3,4),Q(1,4)], [Q(0),Q(0),Q(0),Q(0),Q(1)]]
# Check the split derivative vanishes and the other two forms do not depend on it.
split=[0,0,0,1,-1,0]
assert sum(split[i]*C[0][i][j]*split[j]for i in range(6)for j in range(6))>0
for k in range(3):
 assert all(sum(split[i]*C[k][i][j]*B[j][a]for i in range(6)for j in range(6))==0 for a in range(5))
for k in (1,2):assert all(sum(split[i]*C[k][i][j]for i in range(6))==0 for j in range(6))
Ct=[[[sum(B[i][a]*C[k][i][j]*B[j][b]for i in range(6)for j in range(6))for b in range(5)]for a in range(5)]for k in range(3)]
H=[[u*Ct[0][i][j]+v*Ct[1][i][j]+(1-u-v)*Ct[2][i][j]for j in range(5)]for i in range(5)]
N=u**2*(-20*v**3+420*v**2-224*v-128)+u*(7*v**4+559*v**3-842*v**2+164*v+112)+4*v**4+156*v**3-248*v**2+52*v+36
S=u**2*(156*v**2-84*v-56)+u*(234*v**3-344*v**2+61*v+49)+v**4+65*v**3-102*v**2+20*v+16
M=[[64*x for x in row]for row in H];A=[row+[Poly(1)]for row in M]+[[Poly(1)]*5+[Poly()]]
detM=det(M);adjSum=-det(A)
assert (4*S*detM-64*N*adjSum).iszero()
print('PASS: stationary split and exact T(u,v)=N/(4S) reconstructed from active matrices.',flush=True)
F=N-4*(1-g)*S
fac1=v-1;fac2=104*g*v**2-56*g*v-28*g+v**3-38*v**2+20*v+12
fac3=8328*g*v**3-18144*g*v**2+5628*g*v+3420*g+49*v**4-3367*v**3+7130*v**2-2072*v-1484
Dis=F.coefficient(0,1)**2-4*F.coefficient(0,2)*F.coefficient(0,0)
assert (Dis-fac1*fac2*fac3).iszero()
DP=694*v**3-1512*v**2+469*v+285;GN=-49*v**4+3367*v**3-7130*v**2+2072*v+1484
PT=34006*v**6-148176*v**5+211627*v**4-226450*v**3+422009*v**2-423516*v+105476
assert (GN.deriv(1)*DP-GN*DP.deriv(1)+PT).iszero()
print('PASS: printed discriminant factorization and derivative identity yielding sextic.',flush=True)
# Rigorous interval bounds for the other factors on the displayed isolating box.
class I:
 def __init__(self,l,h=None):self.l=Q(l);self.h=Q(l if h is None else h)
 def __add__(self,o):
  o=o if isinstance(o,I)else I(o);return I(self.l+o.l,self.h+o.h)
 __radd__=__add__
 def __mul__(self,o):
  o=o if isinstance(o,I)else I(o);w=[self.l*o.l,self.l*o.h,self.h*o.l,self.h*o.h];return I(min(w),max(w))
 __rmul__=__mul__
 def __pow__(self,n):return product([self]*n)
vals=[I(0),I('0.359275945611134729','0.359275945611134730'),I('0.41977282048624007772','0.41977282048624007773')]
for f in [fac1,fac2]:
 z=f.ev(vals);assert z.l>0 or z.h<0
print('PASS: the other two discriminant factors are nonzero throughout the isolating box.')