"""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))

# Coefficients of endpoint threshold masses and marginal Q-values.
B=[[[[sum((P[r][k][h]for k,e in enumerate(edges)if i in e and L[e[0]][e[1]]<=s),F(0))for h in range(6)]for s in range(5)]for i in range(6)]for r in range(3)]
S=[[[[sum((B[r][j][L[i][j]][h]for j in range(6)if s<=L[i][j]),F(0))for h in range(6)]for s in range(5)]for i in range(6)]for r in range(3)]
G=[[[2*x for x in row]for row in C]for C in Cs]
strict=0;symbolic=0;min_gap=None;neutrals=[]
for profile in product(range(1,5),repeat=7):
 row=profile[:6];ell=profile[6]
 SG=[[[sum((B[r][j][row[j]][h]for j in range(6)if s<=row[j]),F(0))for h in range(6)]for s in range(5)]for r in range(3)]
 coeff=[[[2*S[r][j][row[j]][h]+2*SG[r][row[j]][h]-G[r][j][h]for h in range(6)]for j in range(6)]+[[2*SG[r][ell][h]for h in range(6)]]for r in range(3)]
 vals=[[sum(c*x for c,x in zip(v,x0[:6]))for v in rr]for rr in coeff];choices=[max(range(7),key=lambda j:vals[r][j])for r in range(3)]
 lower=sum(lam[r]*sum(c*x for c,x in zip(coeff[r][choices[r]],a))for r in range(3))-2*T
 if lower.l>F(19,1000):
  strict+=1;min_gap=lower.l if min_gap is None else min(min_gap,lower.l);continue
 # An exact coefficient identity gives K_r>=gradient_r[oldblock]
 # for each r. Stationarity then gives weighted K>=2T exactly.
 identities=[j for j in range(6)if all(G[r][j]in coeff[r]for r in range(3))]
 assert identities,(profile,float(lower.l))
 symbolic+=1;neutrals.append([list(profile),identities])
assert(strict,symbolic)==(16354,30)
print('Profiles:',strict,'strictly positive,',symbolic,'exactly nonnegative by identities')
print('Uniform strict derivative lower bound >',float(min_gap))
print('Certified gamma stationary interval:',float(1-T.h),float(1-T.l))
print('PASS: no first-order improving added-block profile')
# Neutral profiles and their mutual labels: exact checks for the stronger
# directional, all-orders statement proved in the accompanying note.
neutral_data=[]
for profile,identities in neutrals:
 assert len(identities)==1
 clone=identities[0];row=profile[:6]
 assert all(row[j]>=L[clone][j]for j in range(6))
 assert profile[6]>=L[clone][clone]
 for r in range(3):
  for k,edge in enumerate(edges):
   if clone in edge and any(P[r][k]):
    other=edge[1]if edge[0]==clone else edge[0]
    assert row[other]==L[clone][other]
 SG=[[[sum((B[r][j][row[j]][h]for j in range(6)if s<=row[j]),F(0))for h in range(6)]for s in range(5)]for r in range(3)]
 neutral_data.append((profile,clone,SG))
count=0;min_pair_gain=None
for u,(_,i,sg)in enumerate(neutral_data):
 for _,j,sh in neutral_data[u:]:
  for k in range(1,L[i][j]):
   gains=[sum((2*sg[r][k][h]+2*sh[r][k][h]-G[r][i][h]-G[r][j][h])*a[h]for h in range(6))for r in range(3)]
   gain=max(v.l for v in gains);assert gain>F(9,100)
   count+=1;min_pair_gain=gain if min_pair_gain is None else min(min_pair_gain,gain)
assert count==81
print('Neutral clone domination and all active-flow label agreements: PASS')
print('All',count,'subclone pair labels have strict gain >',float(min_pair_gain))
print('PASS: certificate ingredients for all-orders directional insertion stability')