"""Rational outward-rounded Hopf coefficient enclosure, bound to supplied jets.

Arithmetic approach adapted from the paper workspace's amplitude_enclosure.py,
SHA256 56ec39536b13bc53ad03298c30a2bda43a6213d12e54d2add69e7d39c8460e50.
No floating point enters acceptance. This is not a Lean proof.
"""
from fractions import Fraction as Q
from itertools import product

GRID=2**100
def down(x):return Q((x.numerator*GRID)//x.denominator,GRID)
def up(x):return -down(-x)


class Interval:
    def __init__(self,lower=0,upper=None):
        self.lo=Q(lower);self.hi=Q(lower if upper is None else upper)
        if self.lo>self.hi:raise ValueError('Reversed interval.')
    @staticmethod
    def of(x):return x if isinstance(x,Interval) else Interval(x)
    def __add__(self,other):
        b=self.of(other);return Interval(down(self.lo+b.lo),up(self.hi+b.hi))
    __radd__=__add__
    def __neg__(self):return Interval(-self.hi,-self.lo)
    def __sub__(self,other):return self+-self.of(other)
    def __rsub__(self,other):return self.of(other)+-self
    def __mul__(self,other):
        b=self.of(other);values=[x*y for x in (self.lo,self.hi) for y in (b.lo,b.hi)]
        return Interval(down(min(values)),up(max(values)))
    __rmul__=__mul__
    def reciprocal(self):
        if not (self.lo>0 or self.hi<0):raise ArithmeticError('Interval division includes zero.')
        return Interval(down(1/self.hi),up(1/self.lo))
    def __truediv__(self,other):return self*self.of(other).reciprocal()
    def __rtruediv__(self,other):return self.of(other)*self.reciprocal()
    def record(self):return {'lower':str(self.lo),'upper':str(self.hi)}


class ComplexInterval:
    def __init__(self,real=0,imag=0):self.real=Interval.of(real);self.imag=Interval.of(imag)
    @staticmethod
    def of(x):return x if isinstance(x,ComplexInterval) else ComplexInterval(x)
    def __add__(self,other):
        b=self.of(other);return ComplexInterval(self.real+b.real,self.imag+b.imag)
    __radd__=__add__
    def __neg__(self):return ComplexInterval(-self.real,-self.imag)
    def __sub__(self,other):return self+-self.of(other)
    def __rsub__(self,other):return self.of(other)+-self
    def __mul__(self,other):
        b=self.of(other);return ComplexInterval(self.real*b.real-self.imag*b.imag,self.real*b.imag+self.imag*b.real)
    __rmul__=__mul__
    def conjugate(self):return ComplexInterval(self.real,-self.imag)
    def norm_squared(self):return self.real*self.real+self.imag*self.imag
    def reciprocal(self):
        norm=self.norm_squared()
        if norm.lo<=0:raise ArithmeticError('Complex denominator not separated from zero.')
        return ComplexInterval(self.real/norm,-self.imag/norm)
    def __truediv__(self,other):return self*self.of(other).reciprocal()


class CheckedLinearSolver:
    def __init__(self):self.pivots=[]
    def solve(self,A,b):
        n=len(b)
        if len(A)!=n or any(len(row)!=n for row in A):raise ValueError('Square linear system required.')
        M=[[ComplexInterval.of(z) for z in row]+[ComplexInterval.of(rhs)] for row,rhs in zip(A,b)]
        for k in range(n):
            pivot_row=max(range(k,n),key=lambda i:(M[i][k].real.lo+M[i][k].real.hi)**2+(M[i][k].imag.lo+M[i][k].imag.hi)**2)
            M[k],M[pivot_row]=M[pivot_row],M[k];pivot=M[k][k]
            if pivot.norm_squared().lo<=0:raise ArithmeticError('Uncertified pivot.')
            self.pivots.append({'dimension':n,'column':k,'real':pivot.real.record(),'imaginary':pivot.imag.record(),'modulus_squared_lower':str(pivot.norm_squared().lo)})
            M[k]=[z/pivot for z in M[k]]
            for i in range(n):
                if i==k:continue
                factor=M[i][k];M[i]=[a-factor*b for a,b in zip(M[i],M[k])]
        return [row[-1] for row in M]


def contract(tensor,*vectors):
    n=len(vectors[0]);result=[ComplexInterval() for _ in range(n)]
    for indices,coefficients in tensor.items():
        factor=ComplexInterval(1)
        for vector,i in zip(vectors,indices):factor*=vector[i]
        for i in range(n):result[i]+=factor*ComplexInterval(coefficients[i])
    return result


def certify(model,crossing):
    """The model supplies literal stoichiometry, flux and falling-factorial jets."""
    crossing.exact_certificate(model)
    lo,hi=crossing.isolate(110);t=Interval(lo,hi)
    a1,a2,a3,a4=crossing.coefficients(t);omega_squared=a3/a1
    wl,wh=Q(0),Q(1)
    for _ in range(110):
        mid=(wl+wh)/2
        if mid*mid<omega_squared.lo:wl=mid
        elif mid*mid>omega_squared.hi:wh=mid
        else:break
    if not (wl>0 and wl*wl<=omega_squared.lo and wh*wh>=omega_squared.hi):raise ArithmeticError('Frequency enclosure failed.')
    omega=Interval(wl,wh);inverse_state=[Interval(Q(1,2)),Interval(Q(1,600)),Interval(Q(1,1500)),t/8]
    tensors={order:model.jet(order,inverse_state) for order in (1,2,3)}
    A=[[ComplexInterval(tensors[1][(j,)][i]) for j in range(4)] for i in range(4)]
    M=[[A[i][j]-ComplexInterval(0,omega)*(i==j) for j in range(4)] for i in range(4)];solver=CheckedLinearSolver()
    right=solver.solve([row[:3] for row in M[:3]],[-M[i][3] for i in range(3)])+[ComplexInterval(1)]
    left=solver.solve([[M[j][i] for j in range(3)] for i in range(3)],[-M[3][i] for i in range(3)])+[ComplexInterval(1)]
    pairing=sum((a*b for a,b in zip(left,right)),ComplexInterval());left=[z/pairing for z in left]
    conjugate=[z.conjugate() for z in right];B=tensors[2];C=tensors[3]
    v11=solver.solve(A,contract(B,right,conjugate));v20=solver.solve([[ComplexInterval(0,2*omega)*(i==j)-A[i][j] for j in range(4)] for i in range(4)],contract(B,right,right))
    parts=[contract(C,right,right,conjugate),contract(B,right,v11),contract(B,conjugate,v20)]
    coefficient=sum((left[i]*(parts[0][i]-2*parts[1][i]+parts[2][i]) for i in range(4)),ComplexInterval()).real/(2*omega)
    norm=sum((z.norm_squared() for z in right),Interval());unit_coefficient=coefficient/norm
    if not Q(-23,1000)<unit_coefficient.lo<=unit_coefficient.hi<Q(-22,1000):raise ArithmeticError('Lyapunov sign target not certified.')
    return {'status':'PASS exact rational outward-rounded enclosure; not Lean','bits':100,'t':t.record(),'omega':omega.record(),
        'l1_q4_equal_one':coefficient.record(),'q_norm_squared':norm.record(),'l1_unit_norm':unit_coefficient.record(),'checked_pivots':solver.pivots,
        'binding':'Jets are constructed from the supplied literal S,Y,flux. The exact quartic crossing identity implies det(A-i*omega I)=0; checked leading-block pivots make its Schur complement zero, completing the omitted eigenvector equations. Checked pairing and resolvents justify each division.',
        'scope':'Together with the exact crossing and conventional nondegenerate Hopf theorem this certifies supercritical local attraction. No explicit decimal orbit or size of its parameter neighborhood is validated.'}