"""Resident components adapted from entry 21, with entry-44 readiness and readout."""
from fractions import Fraction as Q
from dataclasses import dataclass, replace
A=Q(1,512000000)

MATRICES = {
    'L': ((1115801,-613708,2346047,3444885),(-613708,1111499,-2339313,-3434384),
          (2346047,-2339313,6767757,10178806),(3444885,-3434384,10178806,15517433)),
    'H': ((846084,-1159490,2352853,3510058),(-1159490,4333988,-6781643,-10233028),
          (2352853,-6781643,11398763,17222276),(3510058,-10233028,17222276,26082111)),
}


@dataclass(frozen=True)
class Interval:
    lo: Q
    hi: Q | None = None
    def __post_init__(self):
        object.__setattr__(self,'lo',Q(self.lo))
        object.__setattr__(self,'hi',Q(self.lo if self.hi is None else self.hi))
        if self.lo>self.hi:raise ValueError('Inverted interval')
    @staticmethod
    def cast(v):return v if isinstance(v,Interval) else Interval(v)
    def __add__(self,v):v=self.cast(v);return Interval(self.lo+v.lo,self.hi+v.hi)
    __radd__=__add__
    def __neg__(self):return Interval(-self.hi,-self.lo)
    def __sub__(self,v):return self+-self.cast(v)
    def __rsub__(self,v):return self.cast(v)+-self
    def __mul__(self,v):
        v=self.cast(v);values=[x*y for x in (self.lo,self.hi) for y in (v.lo,v.hi)]
        return Interval(min(values),max(values))
    __rmul__=__mul__
    def __truediv__(self,v):
        v=self.cast(v)
        if v.lo<=0<=v.hi:raise ZeroDivisionError('Interval denominator includes zero')
        return self*Interval(1/v.hi,1/v.lo)
    def __rtruediv__(self,v):return self.cast(v)/self
    def json(self):return [str(self.lo),str(self.hi)]


def reconstruction(z):
    b=60/(z+2);k=(20004*z*z-159984*z)/20001
    return (z*b+k,b,z,(16*z+2*z*z)/Q(20001,10000))


def residual(z):
    aa,bb,_,_=reconstruction(z)
    return Q(1,100000)*aa*aa+aa+Q(99999,100000)*bb-33


class ChemicalRegions:
    """Rational energy enclosures around tightly isolated stationary roots."""
    def __init__(self,refinements=150):
        self.roots={};self.centers={}
        for tag,left,right in [('L',Q('0.99579401232'),Q('0.99579401233')),
                               ('H',Q('2.97636724376'),Q('2.97636724377'))]:
            fl=residual(left);assert fl*residual(right)<0
            for _ in range(refinements):
                mid=(left+right)/2;fm=residual(mid)
                if fm==0:left=right=mid;break
                if fl*fm<0:right=mid
                else:left=mid;fl=fm
            self.roots[tag]=Interval(left,right)
            self.centers[tag]=reconstruction(self.roots[tag])

    def energy(self,cell):
        delta=[Q(n,cell.size)-center for n,center in zip(cell.counts,self.centers[cell.tag])]
        value=sum((Q(MATRICES[cell.tag][i][j],10**6)*delta[i]*delta[j]
                   for i in range(4) for j in range(4)),Interval(0))
        return Interval(max(Q(0),value.lo),max(Q(0),value.hi))

    def admits(self,cell,threshold,closed=False):
        energy=self.energy(cell)
        if (energy.hi<=threshold if closed else energy.hi<threshold):return True
        if (energy.lo>threshold if closed else energy.lo>=threshold):return False
        raise ArithmeticError('Energy comparison unresolved; refine the stationary-root brackets.')

    def newborn(self,tag,N):
        centers=self.centers[tag]
        def floor_value(v):
            left=(N*v.lo).numerator//(N*v.lo).denominator
            right=(N*v.hi).numerator//(N*v.hi).denominator
            if left!=right:raise ArithmeticError('Stationary floor unresolved; refine root brackets.')
            return left
        cell=Cell(tag,N,tuple(floor_value(v) for v in centers))
        if not self.admits(cell,A,closed=True):
            raise ValueError('Rounded preparation is outside the ready region at this N.')
        return cell


@dataclass(frozen=True)
class Cell:
    tag: str
    size: int
    counts: tuple[int,int,int,int]
    def __post_init__(self):
        if self.tag not in ('H','L') or type(self.size) is not int or self.size<1:
            raise ValueError('A cell needs a valid tag and positive integer size.')
        if len(self.counts)!=4 or any(type(n) is not int or n<0 for n in self.counts):
            raise ValueError('Resident counts must be four nonnegative Python integers.')
    @property
    def readout(self):
        value=self.counts[2]-2*self.size
        return 'H' if value>0 else 'L'


@dataclass(frozen=True)
class Channel:
    name: str
    inputs: tuple[int,...]
    outputs: tuple[int,...]
    coefficient: Q
    def rate(self,cell):
        value=self.coefficient*Q(cell.size)**(1-sum(self.inputs))
        for count,order in zip(cell.counts,self.inputs):
            if count<order:return Q(0)
            for offset in range(order):value*=count-offset
        return value
    def apply(self,cell):
        if self.rate(cell)<=0:raise ValueError('A zero-propensity channel cannot fire.')
        return replace(cell,counts=tuple(n-out+new for n,out,new in zip(cell.counts,self.inputs,self.outputs)))


class ResidentChemistry:
    def __init__(self):
        self.channels=(
            Channel('A to B+z',(1,0,0,0),(0,1,1,0),Q(1)),
            Channel('B+z to A',(0,1,1,0),(1,0,0,0),Q(1)),
            Channel('z to H',(0,0,1,0),(0,0,0,1),Q(16)),
            Channel('H to z',(0,0,0,1),(0,0,1,0),Q(1)),
            Channel('H to 2z',(0,0,0,1),(0,0,2,0),Q(1)),
            Channel('2z to H',(0,0,2,0),(0,0,0,1),Q(2)),
            Channel('feed A',(0,0,0,0),(1,0,0,0),Q(6)),
            Channel('remove A',(1,0,0,0),(0,0,0,0),Q(1)),
            Channel('feed B',(0,0,0,0),(0,1,0,0),Q(27)),
            Channel('remove B',(0,1,0,0),(0,0,0,0),Q(1)),
            Channel('B to 2A',(0,1,0,0),(2,0,0,0),Q(1,100000)),
            Channel('2A to B',(2,0,0,0),(0,1,0,0),Q(1,100000)),
            Channel('remove H',(0,0,0,1),(0,0,0,0),Q(1,10000)),
        )

    def generator(self,cell,observable):
        base=observable(cell)
        return sum((ch.rate(cell)*(observable(ch.apply(cell))-base)
                    for ch in self.channels if ch.rate(cell)),Q(0))