"""Exact reference probabilities and one-sided recurring-division reserves."""
from dataclasses import dataclass
from fractions import Fraction as F
from math import ceil
from exact_engine import reserves,log_enclosure,poly_eval_interval


def boundary(x):
    x=F(x)
    if not 0<x<1:raise ValueError('x=exp(-h) must be strictly between zero and one.')
    A=(x**4+x**3+x*x+x+1)*(3*x**5+6*x**4+9*x**3+12*x*x+8*x+4)
    Q=3*x**8+12*x**7+30*x**6+60*x**5+98*x**4+137*x**3+170*x*x+159*x+66
    return (1-x)*Q/(4*A),4*A*x*x*(1-x)**3/210


def pulse_extinction(x,c,order):
    x,c=F(x),F(c)
    if not 0<x<1 or not -F(1,4)<=c<=F(1,4) or order not in ['AB','BA']:raise ValueError('Invalid exact pulse reference.')
    first,last=((2,4),(3,3)) if order=='AB' else ((3,3),(2,4))
    def B(r,s):return (x**(s+1)-x**(r+s))/(r-1)+(x*x-x**(s+1))/(s-1)
    return 1-x*x-B(first[0],last[0])-B(first[1],last[1])+(F(1,4)+c)*(B(2*first[0],2*last[0])+B(2*first[1],2*last[1]))+(F(1,2)-2*c)*B(sum(first),sum(last))


def constant_extinction(x,c,action):
    x,c=F(x),F(c)
    if action not in ['A','B'] or not 0<x<1 or not -F(1,4)<=c<=F(1,4):raise ValueError('Invalid constant-action reference.')
    r,s=(2,4) if action=='A' else (3,3)
    I=lambda rate:(x-x**rate)/(rate-1)
    return 1-x-I(r)-I(s)+(F(1,4)+c)*(I(2*r)+I(2*s))+(F(1,2)-2*c)*I(r+s)


@dataclass(frozen=True)
class DecisionBand:
    threshold:F
    slope:F
    low_reserve:F
    high_reserve:F
    low_action:str
    high_action:str
    def __post_init__(self):
        if self.slope<=0 or min(self.low_reserve,self.high_reserve)<0:raise ValueError('Positive slope and nonnegative error reserves required.')
    @property
    def low_boundary(self):return self.threshold-self.low_reserve/self.slope
    @property
    def high_boundary(self):return self.threshold+self.high_reserve/self.slope
    def decide(self,c_low,c_high=None):
        c_low=F(c_low);c_high=c_low if c_high is None else F(c_high)
        if not -F(1,4)<=c_low<=c_high<=F(1,4):raise ValueError('Feasible covariance interval required.')
        if c_high<self.low_boundary:return self.low_action
        if c_low>self.high_boundary:return self.high_action
        return 'unresolved'
    def low_minus_high_extinction_interval(self,c_low,c_high=None):
        c_low=F(c_low);c_high=c_low if c_high is None else F(c_high)
        return (self.slope*(self.threshold-c_high)-self.low_reserve,self.slope*(self.threshold-c_low)+self.high_reserve)


class ReferenceCertificate:
    """Fixed response law A=(2,4), B=(3,3); founder division rate one.

    The arbitrary-descendant certificate requires 2*x*x-x**6 >= 1.
    Otherwise only the common half-marginal kernel is covered by this helper.
    """
    def __init__(self,x=F(7,8),epsilon_max=F(1,100),arbitrary_descendants=False):
        self.x=F(x);self.epsilon=F(epsilon_max)
        if not F(1,2)<self.x<1 or self.epsilon<0:raise ValueError('This alternating-log implementation requires 1/2 < exp(-h) < 1 and epsilon >= 0.')
        if arbitrary_descendants and 2*self.x**2-self.x**6<1:raise ValueError('Arbitrary descendant kernels fail the sufficient supersolution criterion.')
        self.arbitrary_descendants=arbitrary_descendants
        self.h=log_enclosure(self.x.denominator,self.x.numerator)
        self.polynomials={name:reserves(name,self.x) for name in ['A','B','AB','BA']}
        self.W={name:poly_eval_interval(p,*self.h) for name,p in self.polynomials.items()}
        self.V={name:reserves(name,self.x,sharp=False)[0] for name in self.polynomials}
        if any(lo<0 or hi>self.V[name] for name,(lo,hi) in self.W.items()):raise ValueError('Failed reserve enclosure.')
    def band(self,task):
        if task=='pulses':
            threshold,slope=boundary(self.x);low,high='AB','BA'
        elif task=='constant':
            a=constant_extinction(self.x,0,'A');b=constant_extinction(self.x,0,'B');slope=4*(constant_extinction(self.x,F(1,4),'A')-a);threshold=(b-a)/slope;low,high='B','A'
        else:raise ValueError('Choose pulses or constant.')
        return DecisionBand(threshold,slope,self.epsilon*self.W[low][1],self.epsilon*self.W[high][1],low,high)
    def reference_gap(self,c,task='pulses'):
        band=self.band(task);return band.slope*(band.threshold-F(c))
    def reversal_range(self,c):
        band=self.band('pulses');c=F(c)
        if not 0<c<=F(1,4):raise ValueError('Positive feasible witness required.')
        low=band.slope*(band.threshold+c);high=band.slope*(c-band.threshold)
        return max(F(0),min(low/self.W['AB'][1],high/self.W['BA'][1]))
    def assay_band(self,task,grid=10**12):
        # Outward coarsening keeps exact binomial denominators manageable.
        b=self.band(task)
        return DecisionBand(b.threshold,b.slope,F(ceil(b.low_reserve*grid),grid),F(ceil(b.high_reserve*grid),grid),b.low_action,b.high_action)