"""Exact local stationarity certificate; standard-library arithmetic only."""
from fractions import Fraction as F
from itertools import combinations,product
from pathlib import Path
import json
HERE=Path(__file__).parent
SOURCE=HERE/'stationarity-input.json'
d=json.loads(SOURCE.read_text());L=d['labels'];Cs=[[[F(x,4)for x in row]for row in C]for C in d['C4']];P=[[[F(x,2)for x in row]for row in p]for p in d['P2']]
x0=list(map(F,d['center']));rad=F(d['radius']);n=9
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,o):o=asI(o);return I(self.l+o.l,self.h+o.h)
 __radd__=__add__
 def __neg__(self):return I(-self.h,-self.l)
 def __sub__(self,o):return self+-asI(o)
 def __rsub__(self,o):return asI(o)+-self
 def __mul__(self,o):
  o=asI(o);v=[self.l*o.l,self.l*o.h,self.h*o.l,self.h*o.h];return I(min(v),max(v))
 __rmul__=__mul__
 def __truediv__(self,o):
  o=asI(o);assert o.l*o.h>0;return self*I(1/o.h,1/o.l)
 def absmax(self):return max(abs(self.l),abs(self.h))
def asI(x):return x if isinstance(x,I)else I(x)
def matvec(A,x):return [sum(c*y for c,y in zip(row,x))for row in A]
def mixed(lam):return [[sum(lam[k]*Cs[k][i][j]for k in range(3))for j in range(6)]for i in range(6)]
def equations(x):
 a=x[:6];lam=[x[6],x[7],1-x[6]-x[7]];H=mixed(lam)
 out=[v-x[8]for v in matvec(H,a)]
 for k in(1,2):out.append(sum(a[i]*(Cs[k][i][j]-Cs[0][i][j])*a[j]for i in range(6)for j in range(6)))
 return out+[sum(a)-1]
def jacobian(x):
 a=x[:6];lam=[x[6],x[7],1-x[6]-x[7]];H=mixed(lam);J=[[0 for _ in range(9)]for _ in range(9)]
 for i in range(6):
  J[i][:6]=H[i];J[i][6]=sum((Cs[0][i][j]-Cs[2][i][j])*a[j]for j in range(6));J[i][7]=sum((Cs[1][i][j]-Cs[2][i][j])*a[j]for j in range(6));J[i][8]=-1
 for k in(1,2):
  for i in range(6):J[5+k][i]=2*sum((Cs[k][i][j]-Cs[0][i][j])*a[j]for j in range(6))
 J[8][:6]=[1]*6;return J
def inverse(A):
 z=[[F(v)for v in row]+[F(i==j)for j in range(len(A))]for i,row in enumerate(A)]
 for k in range(len(A)):
  p=next(i for i in range(k,len(A))if z[i][k]);z[k],z[p]=z[p],z[k];v=z[k][k];z[k]=[x/v for x in z[k]]
  for i in range(len(A)):
   if i!=k:
    v=z[i][k];z[i]=[x-v*y for x,y in zip(z[i],z[k])]
 return [row[len(A):]for row in z]
# Rational Newton preconditioner and interval contraction proof.
R=inverse(jacobian(x0));box=[I(x-rad,x+rad)for x in x0];JI=jacobian(box)
K=[[I(i==j)-sum(R[i][k]*JI[k][j]for k in range(n))for j in range(n)]for i in range(n)]
row_sums=[sum(v.absmax()for v in row)for row in K];rho=max(row_sums)
res=matvec(R,equations(x0));assert rho<1
assert all(abs(res[i])+row_sums[i]*rad<rad for i in range(n))
print('Unique stationary root in certified box; contraction <',float(rho))
a=box[:6];lam=[box[6],box[7],1-box[6]-box[7]];T=box[8]
assert all(x.l>0 for x in a+lam)
# Positive definiteness of the mixed quadratic throughout the box.
H=mixed(lam);LD=[[I(0)for _ in range(6)]for _ in range(6)];diag=[]
for i in range(6):
 pivot=H[i][i]-sum(LD[i][k]*LD[i][k]*diag[k]for k in range(i));assert pivot.l>0;diag.append(pivot);LD[i][i]=I(1)
 for j in range(i+1,6):LD[j][i]=(H[j][i]-sum(LD[j][k]*LD[i][k]*diag[k]for k in range(i)))/pivot
# Each active flow uses six strictly-positive basic coordinates in this box.
edges=list(combinations(range(6),2))+[(i,i)for i in range(6)]
for p in P:
 positive=[i for i,row in enumerate(p)if any(row)]
 assert len(positive)==6
 assert all(matvec(p,a)[i].l>0 for i in positive)
 assert all(sum(p[k][j]for k,e in enumerate(edges)if i in e)==F(i==j)for i in range(6)for j in range(6))
 # Independently reconstruct the quadratic of this matching flow.
 Hflow=[[sum(L[e[0]][e[1]]<=L[i][j] and L[f[0]][f[1]]<=L[i][j]for i in set(e)for j in set(f))for f in edges]for e in edges]
 rec=[[sum(p[u][i]*Hflow[u][v]*p[v][j]for u in range(21)for v in range(21))for j in range(6)]for i in range(6)]
 assert rec==Cs[P.index(p)]
# All-Gram check at the certified stationary root. Enumerate all bases,
# prove the feasibility chamber is constant throughout the root box,
# and compare each symmetrized Gram quadratic to T.
red_edges=[(i,j)for i,j in combinations(range(6),2)if L[i][j]<max(L[i][i],L[j][j])]+[(i,i)for i in range(6)]
m=len(red_edges);AA=[[int(i in e)for e in red_edges]for i in range(6)]
HH=[[sum(L[e[0]][e[1]]<=L[i][j] and L[f[0]][f[1]]<=L[i][j]for i in set(e)for j in set(f))for f in red_edges]for e in red_edges]
pp=[];nb=0
for js in combinations(range(m),6):
 try:inv=inverse([[row[j]for j in js]for row in AA])
 except StopIteration:continue
 nb+=1;vals=matvec(inv,a)
 if any(v.h<0 for v in vals):continue
 assert all(v.l>0 for v in vals),'Root box intersects a basis feasibility wall'
 p2=[[0]*6 for _ in range(m)]
 for j,row in zip(js,inv):
  assert all((2*v).denominator==1 for v in row)
  p2[j]=[int(2*v)for v in row]
 pp.append(p2)
assert nb==2880 and len(pp)==240
assert len({tuple(v for row in p for v in row)for p in pp})==240
hp=[[[sum(HH[u][v]*p[v][j]for v in range(m))for j in range(6)]for u in range(m)]for p in pp]
active_c8=[[[int(8*v)for v in row]for row in C]for C in Cs]
from math import lcm
den=lcm(*(v.denominator for v in x0[:6]));an=[int(v*den)for v in x0[:6]]
num_strict=0;num_active=0;gramgap=None
for s in range(len(pp)):
 for t in range(s,len(pp)):
  mat=[[sum(pp[s][u][i]*hp[t][u][j]for u in range(m))for j in range(6)]for i in range(6)]
  c8=[[mat[i][j]+mat[j][i]for j in range(6)]for i in range(6)]
  if c8 in active_c8:num_active+=1;continue
  q0=F(sum(c8[i][j]*an[i]*an[j]for i in range(6)for j in range(6)),8*den*den)
  # |Q(a)-Q(a0)| <= 2r sum_ij |Cij| |a0j| + r^2 sum_ij |Cij|.
  err=rad*F(sum(abs(c8[i][j])*an[j]for i in range(6)for j in range(6)),4*den)+rad*rad*F(sum(abs(v)for row in c8 for v in row),8)
  gap=T.l-q0-err;assert gap>0
  num_strict+=1;gramgap=gap if gramgap is None else min(gramgap,gap)
assert num_strict+num_active==240*241//2
print('All-Gram root check:',num_strict,'strict pairs;',num_active,'active identities; gap >',float(gramgap))


print("PASS: all matching flows have noncritical density at most the stationary value.")