"""Exact finite-count design decisions and explicit approximation boundaries."""
from fractions import Fraction as F
from model import probability,positive_integer,gate_bound


def batch_certificate(H,m,T,alpha):
    H,alpha=probability(H),probability(alpha);positive_integer(m);positive_integer(T)
    c=(1-H)**T;floor=1-c
    if c<=F(m,m+T):
        power=floor**T;value=float(floor)
    else:
        # Raising to positive integer T removes the irrational exponent m/T.
        power=F(T,m+T)**T*(F(m,m+T)/c)**m
        import math
        value=math.exp((T*math.log(T/(m+T))+m*(math.log(m/(m+T))-math.log(float(c))))/T)
    return dict(floor=floor,bound_numeric=value,certified=power<=alpha**T,
                comparison='Exact comparison of bound**T against alpha**T',reachable=floor<=alpha)


def observation_bound(H,m,ec,et):
    H,ec,et=map(probability,(H,ec,et));positive_integer(m)
    if H==1:return F(1)
    c=1-H
    points=[F(0),F(1),et/c,1-ec,F(m,m+1)*(1+et)/c-ec/(m+1)]
    return max(min(F(1),w+ec)**m*min(F(1),1-c*w+et) for w in points if 0<=w<=1)


def specificity_bound(H,m,kappa):
    H,kappa=map(probability,(H,kappa));positive_integer(m)
    if H==1 or kappa==1:return F(1)
    A=(1-kappa*H)/(1-kappa);Hp=H*(1-kappa)/(1-kappa*H)
    v,z=gate_bound(Hp,m)
    return kappa**m if z<kappa else A*v


def exp_negative(x,degree=160):
    x=F(x)
    if x<0 or type(degree)!=int or degree<2 or degree%2:raise ValueError('Nonnegative argument and positive even degree required.')
    # Range reduction keeps the enclosure informative for large loading values.
    folds=0
    while x>1:x/=2;folds+=1
    term=total=F(1)
    for n in range(1,degree+1):term*=-x/n;total+=term
    lo=max(F(0),total-term*x/(degree+1));hi=min(F(1),total)
    # Round OUTWARD to avoid exponential growth of exact denominators.
    scale=10**40
    def outward(l,u):return F((l*scale).__floor__(),scale),F((u*scale).__ceil__(),scale)
    lo,hi=outward(lo,hi)
    for _ in range(folds):lo,hi=outward(lo*lo,hi*hi)
    return lo,hi


class PoissonReference:
    def __init__(self,loading,a,K):
        self.loading=F(loading);self.a=probability(a);self.K=positive_integer(K)
        if self.loading<=0 or not 0<self.a<=F(1,2):raise ValueError('Positive loading and balanced fraction in (0,.5] required.')

    def chord_status(self):
        lo,hi=exp_negative(2*self.loading);J=(1-2*self.a)**self.K
        if self.K*self.a*(1-hi)>=self.loading*(1-J):return 'certified'
        if self.K*self.a*(1-lo)<self.loading*(1-J):return 'invalid'
        return 'unresolved'

    def chord_bound(self,m):
        if self.chord_status()!='certified':return None
        lo,hi=exp_negative(2*self.loading);g=gate_bound((1-2*self.a)**self.K,m)[0]
        return (1-hi)**m*g,(1-lo)**m*g

    def source_error(self,states,m):
        # states: exact weights and S=X+Y, not independently reloaded per target.
        if sum(F(w) for w,s in states)!=1:raise ValueError('Weights must sum to one.')
        zl=zu=r=F(0)
        for w,s in states:
            w=probability(w);s=F(s)
            if not 0<=s<=2:raise ValueError('S outside [0,2].')
            lo,hi=exp_negative(self.loading*s);zl+=w*(1-hi);zu+=w*(1-lo)
            r+=w*(1-self.a*s)**self.K
        return zl**positive_integer(m)*r,zu**m*r

    def envelope_numeric(self,m):
        """General concave-envelope solution; numerical, never a design certificate.

        A maximizing envelope point corresponds to one or two diagonal source
        states X=Y=S/2. The right affine piece is a tangent to the original curve.
        """
        import math
        from scipy.optimize import brentq,minimize_scalar
        l,a,K=float(self.loading),float(self.a),self.K;d=-math.expm1(-2*l);J=(1-2*a)**K
        if d==1: return {'status':'numerically-unresolved: endpoint rounded to one'}
        def h(z):return (1+a/l*math.log1p(-z))**K
        def hp(z):return -K*a/l*(1+a/l*math.log1p(-z))**(K-1)/(1-z)
        theta=(K-1)*a/l
        if K*a/l >= (1-J)/d:tau=0.
        elif theta<=1-2*a:tau=d
        else:
            turning=-math.expm1((theta-1)*l/a)
            tau=brentq(lambda z:hp(z)-(J-h(z))/(d-z),0,turning)
        slope=(J-h(tau))/(d-tau) if tau<d else 0.
        def env(z):return h(z) if z<=tau else h(tau)+slope*(z-tau)
        candidates=[0.,tau,d]
        for left,right in [(0,tau),(tau,d)]:
            if right>left:candidates.append(minimize_scalar(lambda z:-z**m*env(z),bounds=(left,right),method='bounded',options={'xatol':1e-14}).x)
        z=max(candidates,key=lambda z:z**m*env(z))
        if z<=tau or tau==d:states=[(1.,-math.log1p(-z)/l)]
        else:
            weight=(d-z)/(d-tau);states=[(weight,-math.log1p(-tau)/l),(1-weight,2.)]
        return dict(status='numerical illustration',bound=z**m*env(z),tangent=tau,gate_probability=z,source_states=states)


def constrained_value(H,J,m,p):
    return (2*p-p*p)**m*((1-p)**2+2*p*(1-p)*H+p*p*J)


def constrained_certificate(H,J,m,target,interval_budget=100000):
    H,J,target=map(probability,(H,J,target));positive_integer(m);positive_integer(interval_budget)
    beta,gamma=2*H-2,1-2*H+J;stack=[(F(0),F(1))];used=0
    while stack:
        if used>=interval_budget:return dict(status='budget-exhausted',intervals=used)
        l,r=stack.pop();used+=1;mid=(l+r)/2
        val=constrained_value(H,J,m,mid)
        if val>target:return dict(status='refuted',witness_p=mid,value=val,intervals=used)
        points=[l,r]
        if gamma<0 and l<=-beta/(2*gamma)<=r:points.append(-beta/(2*gamma))
        upper=(2*r-r*r)**m*max(1+beta*p+gamma*p*p for p in points)
        if upper>target:stack.extend([(l,mid),(mid,r)])
    return dict(status='certified',upper=target,intervals=used)