"""Exact polynomial candidate and supplied-certificate interfaces."""
from dataclasses import dataclass
import sympy as s
R=s.Rational

@dataclass(frozen=True)
class Interval:
    lo:object
    hi:object
    def __post_init__(self):
        object.__setattr__(self,'lo',R(self.lo));object.__setattr__(self,'hi',R(self.hi))
        if self.lo>self.hi:raise ValueError('Reversed interval.')
    def __add__(self,b):
        b=as_interval(b);return Interval(self.lo+b.lo,self.hi+b.hi)
    __radd__=__add__
    def __neg__(self):return Interval(-self.hi,-self.lo)
    def __sub__(self,b):return self+-as_interval(b)
    def __rsub__(self,b):return as_interval(b)+-self
    def __mul__(self,b):
        b=as_interval(b);p=[x*y for x in [self.lo,self.hi] for y in [b.lo,b.hi]];return Interval(min(p),max(p))
    __rmul__=__mul__
    def __truediv__(self,b):
        b=as_interval(b)
        if b.lo<=0<=b.hi:raise ValueError('Division interval contains zero.')
        return self*Interval(1/b.hi,1/b.lo)
    def __rtruediv__(self,b):return as_interval(b)/self
    def __pow__(self,n):
        if type(n)is not int or n<0:raise ValueError('Nonnegative integral exponent required.')
        if n==0:return Interval(1,1)
        p=[self.lo**n,self.hi**n];return Interval(0 if n%2==0 and self.lo<=0<=self.hi else min(p),max(p))
    @property
    def magnitude(self):return max(abs(self.lo),abs(self.hi))
    def strings(self):return [str(self.lo),str(self.hi)]

def as_interval(x):return x if isinstance(x,Interval) else Interval(x,x)

def polynomial_box(expr,variables,box):
    if len(box)!=len(variables):raise ValueError('One interval per variable required.')
    answer=Interval(0,0)
    for powers,c in s.Poly(s.expand(expr),*variables,domain=s.QQ).terms():
        term=as_interval(c)
        for power,b in zip(powers,box):term=term*as_interval(b)**int(power)
        answer=answer+term
    return answer

@dataclass(frozen=True)
class ResidualCertificate:
    variables:tuple
    coordinate:int
    alpha:object
    multiplier:object
    coefficients:tuple
    power:int=1
    def evaluate(self,field,box,residual_bounds,multiplier_floor,disturbance=None):
        n=len(self.variables)
        if len(field)!=n or len(self.coefficients)!=n or len(residual_bounds)!=n or len(box)!=n or type(self.power)is not int or self.power<1 or not 0<=self.coordinate<n:raise ValueError('Certificate dimensions, coordinate or power invalid.')
        eps=tuple(map(R,residual_bounds));h0=R(multiplier_floor);disturbance=tuple(map(s.sympify,disturbance if disturbance is not None else [0]*n))
        if len(disturbance)!=n or min(eps)<0 or h0<=0:raise ValueError('Nonnegative residual bounds, positive requested floor, and matching disturbance required.')
        target=(self.variables[self.coordinate]-s.sympify(self.alpha))**self.power*s.sympify(self.multiplier)-sum(q*f for q,f in zip(self.coefficients,field))
        if s.expand(target)!=0:return dict(status='invalid_identity',identity_defect=str(s.expand(target)))
        enclosure=polynomial_box(self.multiplier,self.variables,box);floor=max(enclosure.lo,-enclosure.hi)
        if floor<h0:return dict(status='unresolved_multiplier_floor',multiplier_interval=enclosure.strings(),requested_floor=str(h0))
        Q=[polynomial_box(q,self.variables,box).magnitude for q in self.coefficients];signed=s.expand(sum(q*d for q,d in zip(self.coefficients,disturbance)));eta=polynomial_box(signed,self.variables,box).magnitude
        bound=((sum(q*e for q,e in zip(Q,eps))+eta)/h0)**R(1,self.power)
        return dict(status='certified_conditional',identity_verified=True,multiplier_interval=enclosure.strings(),floor=str(h0),coefficient_bounds=list(map(str,Q)),signed_load=str(signed),signed_load_bound=str(eta),absolute_error_bound=str(bound),scope='Applies only at states in the supplied box with the declared perturbed-field residual bounds.')

class BlockCandidates:
    """Explicit lexicographic block order: all other species before target.
    Exact positive real roots are candidates, never automatic ACR decisions.
    """
    def __init__(self,variables,target):
        if target not in variables:raise ValueError('Target must be a model species.')
        self.target=target;self.other=tuple(x for x in variables if x!=target)
    def compute(self,field):
        G=s.groebner(field,*self.other,self.target,order='lex',domain=s.QQ);basis=[p.as_expr() for p in G.polys]
        if basis==[1]:return dict(status='empty_algebraic_locus',basis=['1'],candidates=[])
        univariate=[p for p in basis if not p.free_symbols.intersection(self.other)]
        coeffs=univariate or [s.Poly(p,*self.other).terms()[0][1] for p in basis]
        roots=set()
        for c in coeffs:
            for r in s.Poly(c,self.target).real_roots():
                if r>0:roots.add(r)
        return dict(status='candidate_list_only',order=[str(x) for x in self.other]+[str(self.target)],basis=list(map(str,basis)),leading_coefficients=list(map(str,coeffs)),candidates=list(map(str,sorted(roots,key=lambda x:float(x.evalf())))),scope='Completeness for zero-divisor values follows from the paper. A positive regular steady state is additionally sufficient for coverage of every ACR value; neither condition proves that a candidate is ACR.')