#!/usr/bin/env python3
"""DAE certificates: exact Fraction linear algebra/polynomials and 80-digit checks.
Requires mpmath. Numerical checks are diagnostics, not interval enclosures.
No hidden file inputs; --output writes one JSON, otherwise prints it.
Finite tests supplement the proofs; they do not prove arbitrary-DAE assertions.
"""
from fractions import Fraction as F
from collections import Counter
from pathlib import Path
from math import factorial
import argparse,json
import mpmath as mp
mp.mp.dps=80

def eye(n):return [[F(i==j) for j in range(n)] for i in range(n)]
def zeros(n,m):return [[F(0) for _ in range(m)] for _ in range(n)]
def tr(A):return list(map(list,zip(*A)))
def add(A,B):return [[x+y for x,y in zip(a,b)] for a,b in zip(A,B)]
def scale(a,A):return [[a*x for x in row] for row in A]
def mm(A,B):return [[sum(x*y for x,y in zip(a,b)) for b in zip(*B)] for a in A]
def mv(A,v):return [sum(a*x for a,x in zip(row,v)) for row in A]
def power(A,n):
 R=eye(len(A))
 for _ in range(n):R=mm(R,A)
 return R

def rref(A):
 A=[[F(x) for x in row] for row in A];m=len(A);n=len(A[0]);row=0;piv=[]
 for col in range(n):
  pivot=next((i for i in range(row,m) if A[i][col]),None)
  if pivot is None:continue
  A[row],A[pivot]=A[pivot],A[row];a=A[row][col];A[row]=[x/a for x in A[row]]
  for i in range(m):
   if i!=row:
    a=A[i][col];A[i]=[x-a*y for x,y in zip(A[i],A[row])]
  piv.append(col);row+=1
  if row==m:break
 return A,piv

def inv(A):
 n=len(A);R,p=rref([a+b for a,b in zip(A,eye(n))])
 if p[:n]!=list(range(n)):raise ValueError('singular matrix')
 return [row[n:] for row in R]
def det(A):
 A=[[F(x) for x in row] for row in A];d=F(1);n=len(A)
 for j in range(n):
  k=next((k for k in range(j,n) if A[k][j]),None)
  if k is None:return F(0)
  if k!=j:A[k],A[j]=A[j],A[k];d=-d
  pivot=A[j][j];d*=pivot
  for k in range(j+1,n):
   factor=A[k][j]/pivot
   for l in range(j+1,n):A[k][l]-=factor*A[j][l]
 return d

def null_columns(A):
 R,p=rref(A);n=len(A[0]);cols=[]
 for j in range(n):
  if j in p:continue
  v=[F(0)]*n;v[j]=F(1)
  for i,k in enumerate(p):v[k]=-R[i][j]
  cols.append(v)
 return cols

def fitting_certificate(E,A):
 # sigma=0 is valid for the generated examples; input A is explicitly checked.
 d=len(E);B=mm(inv(A),E);Bd=power(B,d);_,piv=rref(Bd)
 image=[list(col) for j,col in enumerate(zip(*Bd)) if j in piv]
 kernel=null_columns(Bd);Q=tr(image+kernel);r=len(image);m=len(kernel)
 T=mm(mm(inv(Q),B),Q);L=zeros(d,d)
 if r:
  B1=[row[:r] for row in T[:r]];Bi=inv(B1)
  for i in range(r):
   for j in range(r):L[i][j]=Bi[i][j]
 for i in range(r,d):L[i][i]=1
 P=mm(mm(L,inv(Q)),inv(A))
 return P,Q,r,m

def poly_derivative(coeff):return [F(k)*coeff[k] for k in range(1,len(coeff))] or [F(0)]
def poly_eval(coeff,t):return sum(a*t**j for j,a in enumerate(coeff))
def mix_polys(A,cs):
 D=max(map(len,cs));return [[sum(A[i][j]*(cs[j][k] if k<len(cs[j]) else 0) for j in range(len(cs))) for k in range(D)] for i in range(len(A))]
def poly_add(a,b):return [(a[k] if k<len(a) else F(0))+(b[k] if k<len(b) else F(0)) for k in range(max(len(a),len(b)))]
def poly_equal(a,b):return all(x==0 for x in poly_add(a,[-x for x in b]))
def mpf(x):return mp.mpf(x.numerator)/x.denominator if isinstance(x,F) else mp.mpf(x)
def fmt(x):return mp.nstr(x,45)

def run():
 c=Counter()
 def ck(name,value):
  if not value:raise ArithmeticError((name,c[name]))
  c[name]+=1
 tol=mp.mpf('1e-65')
 # Coupled algebraic initialization over a rational family, with explicit g_t.
 for u in [F(3,2),F(2),F(5,2),F(3),F(4)]:
  x=u+1/u;z=[u,1/u];xd=-1/u;G=[[F(1),F(1)],[1/u,u]]
  zd=mv(inv(G),[xd+1,0]);ck('index1_initial_constraint',z[0]+z[1]-x==0 and z[0]*z[1]==1)
  ck('index1_initial_derivative',sum(zd)==xd+1 and zd[0]/u+u*zd[1]==0)
  for ih in range(1,21):
   h=mp.mpf(ih)/40;aa=mpf(x)+h+h*h;disc=aa*aa-4*(1+h)
   ck('index1_branch_discriminant',disc>0)
   Z=(aa+mp.sqrt(disc))/2;V=1/Z;X=(mpf(x)+h*Z)/(1+h)
   ck('index1_backward_euler_differential',abs(X-mpf(x)-h*(-X+Z))<tol)
   ck('index1_backward_euler_constraints',abs(Z+V-X-h)<tol and abs(Z*V-1)<tol)
   ck('index1_branch_order',Z>V)
   # Exact determinant formula at rational trial states, independent of roots.
   H=F(ih,40);a=F(ih+11,7);b=F(5,ih+2)
   J=[[1+H,-H,0],[-1,1,1],[0,b,a]]
   ck('index1_newton_determinant',det(J)==a-(1+H)*b)
 ck('terminal_initial_derivative',mv(inv([[1,1],[F(1,2),2]]),[F(1,2),0])==[F(2,3),F(-1,6)])
 s73=mp.sqrt(73);Z1=(13+s73)/8;Z2=(13-s73)/12;X=(53+s73)/24
 ck('terminal_coupled_step',abs(3*X-5-Z1)<tol and abs(Z1+Z2-X-mp.mpf(1)/2)<tol and abs(Z1*Z2-1)<tol)
 ck('terminal_step_determinant',abs(Z1-mp.mpf(3)/2*Z2-s73/4)<tol)
 # Exact quadratic-field certificates for sqrt(73), beyond a decimal root check.
 def qa(a,b):return (a[0]+b[0],a[1]+b[1])
 def qm(a,b):return (a[0]*b[0]+73*a[1]*b[1],a[0]*b[1]+a[1]*b[0])
 def qs(k,a):return (k*a[0],k*a[1])
 ze=(F(13,8),F(1,8));ve=(F(13,12),F(-1,12));xe=(F(53,24),F(1,24))
 ck('terminal_exact_quadratic_product',qm(ze,ve)==(1,0))
 ck('terminal_exact_quadratic_sum',qa(ze,ve)==qa(xe,(F(1,2),F(0))))
 ck('terminal_exact_quadratic_step',qs(3,xe)==qa((F(5),F(0)),ze))
 # At derivative-array order k, vary jets while the state is fixed. The first
 # derivative is determined exactly when all null vectors have zero first block.
 for nu in range(1,6):
  N=zeros(nu,nu)
  for i in range(nu-1):N[i][i+1]=1
  for order in range(nu+1):
   D=zeros(nu*(order+1),nu*(order+1))
   for j in range(order+1):
    for a in range(nu):
     for b in range(nu):D[j*nu+a][j*nu+b]=N[a][b]
     if j:D[j*nu+a][(j-1)*nu+a]=-1
   determined=all(all(x==0 for x in vec[:nu]) for vec in null_columns(D))
   ck('nilpotent_derivative_array_minimality',determined==(order>=nu))
 # General regular pencils, exact polynomial solutions, and a SECOND computed
 # decomposition from kernel/image of (A^-1 E)^d, rather than the supplied P,Q.
 models=0
 for r in [0,1,2]:
  for nu in range(0,6):
   d=r+nu
   if not d:continue
   Ec=zeros(d,d);Ac=eye(d)
   for i in range(r):Ec[i][i]=1;Ac[i][i]=-i-1
   for i in range(r,d-1):Ec[i][i+1]=1
   P=eye(d);Q=eye(d)
   for i in range(1,d):P[i][i-1]=F(i+1,2);Q[i-1][i]=F(2,i+1)
   E=mm(mm(inv(P),Ec),inv(Q));A=mm(mm(inv(P),Ac),inv(Q))
   ck('regular_supplied_left_right',mm(mm(P,E),Q)==Ec and mm(mm(P,A),Q)==Ac)
   ck('regular_invertible_factors',det(P)!=0 and det(Q)!=0 and det(A)!=0)
   # Construct canonical polynomial forcing, not merely a single trigonometric example.
   ys=[[F((-1)**(i+k)*(i+k+1),k+1) for k in range(5)] for i in range(r)]
   a=[poly_add(poly_derivative(ys[i]),[(i+1)*v for v in ys[i]]) for i in range(r)]
   b=[[F((-1)**(i+k)*(i+1),k+2) for k in range(nu+3)] for i in range(nu)]
   z=[]
   for i in range(nu):
    zi=[F(0)]
    for j in range(nu-i):
     q=b[i+j]
     for _ in range(j):q=poly_derivative(q)
     zi=poly_add(zi,[-v for v in q])
    z.append(zi)
   xc=mix_polys(Q,ys+z);fc=mix_polys(inv(P),a+b)
   lhs=mix_polys(E,[poly_derivative(q) for q in xc]);rhs=[poly_add(p,q) for p,q in zip(mix_polys(A,xc),fc)]
   for p,q in zip(lhs,rhs):ck('regular_polynomial_trajectory_identity',poly_equal(p,q))
   PP,QQ,rr,mm0=fitting_certificate(E,A);EE=mm(mm(PP,E),QQ);AA=mm(mm(PP,A),QQ)
   ck('fitting_free_dimension',rr==r and mm0==nu)
   ck('fitting_invertibility',det(PP)!=0 and det(QQ)!=0)
   for i in range(d):
    for j in range(d):
     if i<r and j<r:ck('fitting_differential_identity',EE[i][j]==F(i==j))
     elif i>=r and j>=r:ck('fitting_algebraic_identity',AA[i][j]==F(i==j))
     else:ck('fitting_cross_blocks',EE[i][j]==0 and AA[i][j]==0)
   if nu:
    N=[row[r:] for row in EE[r:]];ck('fitting_nilpotent_exact_index',power(N,nu)==zeros(nu,nu) and power(N,nu-1)!=zeros(nu,nu))
   for sigma in [F(1),F(3,2),F(7)]:
    ck('regular_pencil_determinants',det(add(scale(sigma,E),scale(-1,A)))==det(add(scale(sigma,Ec),scale(-1,Ac)))/(det(P)*det(Q)))
   models+=1
 # Terminal mixed four-dimensional pencil and initial data.
 P=[[1,0,0,0],[1,1,0,0],[0,1,1,0],[0,0,1,1]];Q=[[1,1,0,0],[0,1,1,0],[0,0,1,1],[0,0,0,1]]
 E=[[1,-1,1,-1],[-1,1,0,0],[1,-1,0,1],[-1,1,0,-1]];A=[[-1,1,-1,1],[1,0,0,0],[-1,0,1,-1],[1,0,-1,2]]
 Ec=[[1,0,0,0],[0,0,1,0],[0,0,0,1],[0,0,0,0]];Ac=[[(-1 if i==0 else 1) if i==j else 0 for j in range(4)] for i in range(4)]
 ck('terminal_mixed_certificate',mm(mm(P,E),Q)==Ec and mm(mm(P,A),Q)==Ac and det(P)==det(Q)==1)
 ck('terminal_initial_vector',mv(E,[-8,-6,0,0])==mv(A,[2,0,0,0]))
 for t in [mp.mpf(k)/8 for k in range(-8,25)]:
  fn=lambda u:[-5*u-1+3*mp.exp(-u),-6*u-3*u*u,-3*u*u-u**3,-u**3]
  xx=fn(t);dd=[mp.diff(lambda u:fn(u)[i],t) for i in range(4)];ff=[t,-t,t,t**3-t]
  for a,b in zip(mv(E,dd),[x+y for x,y in zip(mv(A,xx),ff)]):ck('terminal_mixed_trajectory',abs(a-b)<tol)
 for omega in [2,4,8,16,32,64,128]:
  fn=lambda t:[omega*mp.sin(omega*t),omega*mp.sin(omega*t)-mp.cos(omega*t),-mp.cos(omega*t)-mp.sin(omega*t)/omega,-mp.sin(omega*t)/omega]
  ck('high_frequency_consistent_initial_shift',fn(0)==[0,-1,-1,0])
  t=mp.pi/(2*omega);ck('high_frequency_small_data_large_state',0<t<=1 and abs(fn(t)[0]-omega)<tol)
  for t in [mp.mpf(k)/11 for k in range(12)]:
   dd=[mp.diff(lambda u:fn(u)[i],t) for i in range(4)];forcing=[0,0,0,mp.sin(omega*t)/omega]
   for a,b in zip(mv(E,dd),[x+y for x,y in zip(mv(A,fn(t)),forcing)]):ck('high_frequency_original_equations',abs(a-b)<tol)
 # Exact backward Euler initialization impulses and residual amplification.
 for nu in range(1,6):
  N=zeros(nu,nu)
  for i in range(nu-1):N[i][i+1]=1
  for h in [F(1,2),F(1,4),F(1,10),F(2)]:
   K=inv(add(eye(nu),scale(-1/h,N)));finite=zeros(nu,nu)
   for j in range(nu):finite=add(finite,scale(h**(-j),power(N,j)))
   ck('backward_euler_finite_inverse',K==finite)
   prev=[F(0)]*(nu-1)+[F(1)]
   for step in range(nu):
    new=mv(K,mv(scale(-1/h,N),prev))
    ck('backward_euler_impulse_exact_residual',mv(N,[(a-b)/h for a,b in zip(new,prev)])==new);prev=new
   ck('backward_euler_impulse_finite_death',prev==[0]*nu)
   residual=[F(0)]*(nu-1)+[F(1,1000)];error=mv(scale(-1,K),residual)
   ck('backward_euler_residual_to_error',mv(add(scale(1/h,N),scale(-1,eye(nu))),error)==residual)
   ck('backward_euler_chain_residual_gain',abs(error[0])==F(1,1000)*h**(-(nu-1)))
 # Two independent nonlinear constraints and non-diagonal SPD mass matrices.
 cases=0
 for a0 in [F(1,2),F(1),F(3,2)]:
  for b0 in [F(-1),F(1,3),F(2)]:
   q=[a0,b0,F(2,5)];v=[F(3,7),F(-2,5),F(1,11)];G=[q,[0,0,1]]
   L=[[2,1,0],[0,3,1],[1,0,2]];M=add(mm(tr(L),L),eye(3));Mi=inv(M)
   W=mm(mm(G,Mi),tr(G));g=[(sum(x*x for x in q)-1)/2,q[2]-F(1,3)];w=mv(G,v);H=[sum(x*x for x in v),0];force=[F(2),F(-3),F(1,2)]
   ck('multiple_constraints_gram_positive',W[0][0]>0 and det(W)>0)
   for alpha,beta in [(F(1),F(1)),(F(2),F(3)),(F(4),F(1,2))]:
    lam=mv(inv(W),[a+b+2*alpha*c0+beta*beta*d for a,b,c0,d in zip(mv(G,mv(Mi,force)),H,w,g)])
    acc=mv(Mi,[a-b for a,b in zip(force,mv(tr(G),lam))])
    ck('multiple_constraints_dynamic_row',add([[x] for x in mv(M,acc)],[[x] for x in mv(tr(G),lam)])==[[x] for x in force])
    epp=[a+b for a,b in zip(mv(G,acc),H)]
    ck('multiple_constraints_feedback_row',all(a+2*alpha*b+beta*beta*d==0 for a,b,d in zip(epp,w,g)))
    cases+=1
 # The two terminal ellipse states, including mass-weighted force and signs.
 M=[[2,0],[0,3]];Mi=inv(M);force=[F(2),F(-1)];v=[F(-4,5),F(1,5)];H=F(18,25)
 for q,expected_lambda,expected_acc in [([F(1,3),F(2,3)],F(822,875),[F(738,875),F(-657,875)]),([F(2,5),F(4,5)],F(11,4),[F(9,20),F(-9,5)])]:
  G=[[q[0],2*q[1]]];g=(q[0]**2+2*q[1]**2-1)/2;w=mv(G,v)[0];W=mm(mm(G,Mi),tr(G))[0][0]
  lam=(mv(G,mv(Mi,force))[0]+H+6*w+9*g)/W;acc=mv(Mi,[a-b*lam for a,b in zip(force,G[0])])
  ck('terminal_ellipse_multiplier',lam==expected_lambda and acc==expected_acc)
  ck('terminal_ellipse_feedback',mv(G,acc)[0]+H+6*w+9*g==0)
  for h in [F(1,10),F(1,3),F(2,3)]:
   qp=[a+h*b for a,b in zip(q,v)];gp=(qp[0]**2+2*qp[1]**2-1)/2
   ck('ellipse_physical_euler_drift',gp==g+h*w+h*h*H/2)
 partial=sum((F(6)**k)/factorial(k) for k in range(14))
 ck('terminal_exp_lower_certificate',partial==F(10059173,25025) and partial>400)
 rational=F(1,100000)+F(153993,100000)*F(1,400)
 ck('terminal_tolerance_rational_certificate',rational==F(154393,40000000) and rational<F(1,250))
 # Forced critical residual: analytic complex response vs independent quadrature.
 for omega in [mp.mpf(1)/2,mp.mpf(2),mp.mpf(3)]:
  for kappa in [mp.mpf(1)/3,mp.mpf(1),mp.mpf(4)]:
   for t in [mp.mpf(1)/7,mp.mpf(1)/2,mp.mpf(1),mp.mpf(2),mp.mpf(4)]:
    response=lambda u:(mp.exp(1j*kappa*u)-(1+(omega+1j*kappa)*u)*mp.exp(-omega*u))/(omega+1j*kappa)**2
    value=response(t);direct=mp.quad(lambda s:(t-s)*mp.exp(-omega*(t-s))*mp.exp(1j*kappa*s),[0,t/2,t])
    bound=(1-(1+omega*t)*mp.exp(-omega*t))/omega**2
    ck('critical_forcing_convolution_identity',abs(value-direct)<tol)
    ck('critical_forcing_differential_identity',abs(mp.diff(response,t,2)+2*omega*mp.diff(response,t)+omega**2*value-mp.exp(1j*kappa*t))<tol)
    ck('critical_forcing_uniform_bound',abs(value)<=bound+tol)
 # Formal-page constant-defect example, independently differentiated.
 efun=lambda t:mp.mpf(1)/1000+mp.mpf(9)/1000*(1+2*t)*mp.exp(-2*t)
 for t in [mp.mpf(k)/3 for k in range(13)]:
  ck('constant_defect_exact_equation',abs(mp.diff(efun,t,2)+4*mp.diff(efun,t)+4*efun(t)-mp.mpf(1)/250)<tol)
 ck('constant_defect_initial_state',abs(efun(0)-mp.mpf(1)/100)<tol and abs(mp.diff(efun,0))<tol)
 for omega in [F(1,2),F(1),F(3)]:
  B=[[0,1],[-omega**2,-2*omega]];C=add(B,scale(omega,eye(2)));ck('critical_jordan_nilpotency',power(C,2)==zeros(2,2) and C!=zeros(2,2))
  for sh in [F(1,4),F(1,2),F(1),F(3,2),F(2),F(5,2),F(3)]:
   h=sh/omega;r=1-sh;R=add(eye(2),scale(h,B))
   for n in range(1,19):
    predicted=add(scale(r**n,eye(2)),scale(n*h*r**(n-1),C))
    ck('critical_euler_jordan_power',power(R,n)==predicted)
   if sh==1:ck('critical_euler_zero_eigenvalue',power(R,2)==zeros(2,2))
   if sh==2:
    for n in range(1,12):ck('critical_euler_boundary_growth',mv(power(R,n),[1,0])[0]==(-1)**(n-1)*(2*n-1))
 ck('terminal_boundary_fourth_step',F(11,50)*(-7)==F(-77,50))
 ef=mp.mpf(77)/50*mp.exp(-6);bf=mp.mpf(1)/100000+mp.mpf(153993)/100000*mp.exp(-6)
 return {'status':'PASS','checks':sum(c.values()),'categories':dict(c),'precision_digits':80,'regular_pencil_models':models,'multiple_constraint_cases':cases,'terminal':{'backward_euler':{'X':fmt(X),'Z1':fmt(Z1),'Z2':fmt(Z2)},'feedback_exact_e_at_2':fmt(ef),'feedback_defect_bound_at_2':fmt(bf),'rational_bound':str(rational),'tolerance':'1/250'},'scope':'Exact finite certificates and 80-digit diagnostics; continuous existence, index minimality and infinite-time stability rely on the written proofs.'}

if __name__=='__main__':
 p=argparse.ArgumentParser();p.add_argument('--output',type=Path);a=p.parse_args();result=run();s=json.dumps(result,ensure_ascii=False,indent=2)+'\n'
 if a.output:a.output.parent.mkdir(parents=True,exist_ok=True);a.output.write_text(s)
 else:print(s,end='')