"""Fresh scalar certificates for the paper's fixed reference source.

These reproduce finite inequalities, not the general probability theorems or
Lean build. Directed integration follows manuscript Appendix D/check_paper.py.
"""
from fractions import Fraction as F
from math import factorial, prod, floor, ceil, log
from decimal import Decimal as D, localcontext, ROUND_FLOOR, ROUND_CEILING
from scipy.special import digamma, gammaincc, gammainc, gammaln, kv
from scipy.optimize import brentq
import numpy as np
from source import exp_negative


def timing_dual_certificate():
    q=[F(5,500),F(404,500),F(603,500),F(602,500),F(401,500)]
    def c(n,i):
        return F(0) if i<n else prod(q[n:])/q[i]/prod(q[j]-q[i] for j in range(n,5) if j!=i)
    ep=exp_negative(F(4));coef=[]
    for i in [0,4,1,3,2]:
        B=sum(F(4)**n/factorial(n)*c(n,i) for n in range(5))
        vals=[q[i]*(e*B-5*c(0,i)) for e in ep];scale=10**12
        coef.append((F(floor(min(vals)*scale),scale),F(ceil(max(vals)*scale),scale)))
    def params(side):
        z=[c[side] for c in coef];return -z[0],z[1],-z[2],z[3],-z[4]
    def H(B,Dd,E,G,x):return B-Dd*x**3+x**201*(E-G*x)
    A,B,Dd,E,G=params(1);v=F(9901,10000)
    upper=v**396*H(B,Dd,E,G,v)-A
    assert min(Dd,G,396*B-399*Dd*v**3,597*E-598*G*v)>=0 and upper<0
    A,B,Dd,E,G=params(0);grid=list(map(F,['.992','.994','.996','.998','1']));cells=[]
    for l,u in zip(grid,grid[1:]):
        derivative=-3*Dd*l*l+201*u**200*(E-G*l)-G*l**201
        lower=l**396*H(B,Dd,E,G,u)-A
        assert min(Dd,G,E-G*u,H(B,Dd,E,G,u),lower)>0 and derivative<0
        cells.append(dict(left=l,right=u,derivative_upper=derivative,polynomial_lower=lower))
    # x=exp(-t/500): t>=5 gives x<.9901; t<=4 gives x>.992.
    assert exp_negative(F(1,100))[1]<v and exp_negative(F(1,125))[0]>grid[0]
    return dict(coefficients=coef,upper_polynomial_at_v=upper,cells=cells,
                conclusion='loaded hit density >=5 times blank density before 4; <= after 5; manuscript likelihood-ratio theorem connects the middle interval')


def _decimal(x,rounding):
    with localcontext() as ctx:
        ctx.prec=50;ctx.rounding=rounding
        return D(x.numerator)/D(x.denominator)


def _exp(x,upper):
    # Decimal exp is correctly rounded; one next representable endpoint encloses it.
    x=_decimal(x,ROUND_FLOOR if upper else ROUND_CEILING)
    with localcontext() as ctx:
        ctx.prec=50
        y=(-x).exp(context=ctx)
        return y.next_plus(context=ctx) if upper else y.next_minus(context=ctx)


def large_threshold_certificate():
    """Reference a=.01, lambda=4, h>=10^6; no giant generator allocated."""
    with localcontext() as ctx:
        ctx.prec=50;ctx.rounding=ROUND_FLOOR
        l=F(1,20);integral=D(0)
        while l<10:
            u=min(F(10),l+F(1,200))
            integral+=_decimal((u-l)/u,ROUND_FLOOR)*_exp(max(l+F(31,200)/l,u+F(31,200)/u),False)
            l=u
        product_lower=D('.0097')*integral
    with localcontext() as ctx:
        ctx.prec=50;ctx.rounding=ROUND_CEILING
        l=F(293,1000);e1=D(0)
        while l<20:
            u=min(F(20),l+F(1,200));e1+=_decimal((u-l)/l,ROUND_CEILING)*_exp(l,True);l=u
        e1+=_exp(F(20),True)/D(20)
        reserve_blank=D('.01')*e1+_decimal(F(1,9900),ROUND_CEILING)*_exp(F(293,1000),True)
    el,eu=exp_negative(F(4));h=10**6
    dm=4*exp_negative(F(h,60300))[1]+2*exp_negative(F(h,80602))[1]
    dp=2*exp_negative(F(h,30000))[1]+2*exp_negative(F(h,20200))[1]
    loaded_lower=el*sum(F(4)**n/factorial(n)*F(145,1000)/(n+F(155,1000)) for n in range(6))
    z=F(307,250)
    reserve_miss=eu*(sum(z**n/factorial(n)**2 for n in range(21))+z**21/factorial(21)**2/(1-z/F(22**2)))
    poisson_tail=eu*F(4)**26/factorial(26)/(1-F(4,27))
    reciprocal_sum=eu*sum(F(4)**n/factorial(n)/(n-1) for n in range(2,26))+poisson_tail
    logerr=F(301,10000);v=F(129,25000);sqrtv=F(449,6250)
    assert dm<F(1,100000) and dp<F(1,10**12) and poisson_tail<F(3,10**13)
    assert 3*F(1,100)+F(25,h-24)<logerr
    assert F(3,20)/(1-logerr)<F(31,200) and F(3,20)*(1-logerr)>F(29,200)
    assert 2*F(1,100)+F(25,h-24)<F(201,10000)
    assert F(3,10)/(1-F(201,10000))<F(307,1000) and F(3,10)*(1-F(201,10000))>F(293,1000)
    assert eu<F(19,1000) and 4*eu<F(733,10000) and reciprocal_sum<F(39,100)
    assert F(1,200)/(1-logerr)<v and F(1,200)*(1-logerr)>F(1,250) and sqrtv**2>v
    occupied_miss=2*F(733,10000)*sqrtv+v*F(39,100)+F(1,100000)+F(3,10**13)
    assert occupied_miss<F(63,5000)
    # Integral tail/log bound: exp(8.4)>2501, verified without rounded logarithms.
    assert exp_negative(F(42,5))[1]<F(1,2501)
    blank_upper=F(84,1000)+F(1,9900)+F(1,20000)+F(1,100000)
    assert exp_negative(F(8))[1]<F(1,2500) and exp_negative(F(1))[1]<F(2,5) and exp_negative(F(10))[1]<F(1,20000)
    blank_lower=F(product_lower)-dm;miss_lower=loaded_lower-dm
    assert blank_lower>F(1086,100000) and miss_lower>F(563,10000)
    reserve_B=F(reserve_blank)+dp;reserve_M=reserve_miss+dp+poisson_tail
    assert reserve_B<F(9373,10**6) and reserve_M<F(48732,10**6)
    joint_M=1-F(99,100)*(1-F(19,1000)-F(63,5000))
    assert joint_M==F(10321,250000) and blank_upper/100<F(85,100000)
    return dict(minimum_threshold=h,product_integral_lower=str(product_lower),loaded_product_lower=loaded_lower,
                delta_minus=dm,delta_plus=dp,poisson_tail=poisson_tail,
                separating_blank_lower=blank_lower,separating_miss_lower=miss_lower,
                identity_blank_upper=blank_upper/100,identity_miss_upper=joint_M,
                reserve_blank_upper=reserve_B,reserve_miss_upper=reserve_M,
                evidence='Directed scalar bounds; gamma/product comparison and concentration theorems imported from manuscript; Lean not rerun')


def explicit_deadlines(h):
    if type(h)!=int or h<10**6:raise ValueError('Uniform theorem requires integer h>=10^6.')
    a=.01;k=h//2;ell=h-k;c=h/(h+a)
    # Digamma differences evaluate finite reciprocal sums in O(1) storage.
    mu_minus=digamma(h+1)-digamma(h-k+1)+digamma(h+a)-digamma(k+a)
    A=log(k+a)+log(ell+1)+mu_minus
    # z/[(z+a)(2h-z)] = -a/(2h+a)/(z+a) + 2h/(2h+a)/(2h-z).
    mu_plus=-a/(2*h+a)*(digamma(h+a)-digamma(a))+2*h/(2*h+a)*(digamma(2*h+1)-digamma(h+1))
    return dict(threshold=h,zero_headroom_separation=c*(A-log(.15)),identity=c*(A+log(200)),double_capacity=log(h+a)+mu_plus-log(.3),evidence='numerical deadline evaluation, normalized g=1')


class LimitLaw:
    """Numerical infinite-threshold law, not a finite-threshold certificate."""
    def __init__(self,a=.01,loading=4.,headroom_zero=False):
        if not np.isfinite(a) or a<=0 or not np.isfinite(loading) or loading<0:raise ValueError('Invalid limit parameters.')
        self.a,self.loading,self.headroom_zero=a,loading,headroom_zero

    def errors(self,y):
        if not np.isfinite(y) or y<=0:raise ValueError('Positive transformed deadline required.')
        from scipy.stats import poisson
        n=np.arange(60);s=n+self.a;weights=poisson.pmf(n,self.loading)
        if poisson.sf(59,self.loading)>1e-12:raise ValueError('Increase explicit Poisson truncation budget.')
        if self.headroom_zero:
            with np.errstate(over='ignore',invalid='ignore',under='ignore'):
                tail=2*np.exp(s/2*np.log(y)-gammaln(s))*kv(s,2*np.sqrt(y))
            if np.any(~np.isfinite(tail)):
                import mpmath as mp
                with mp.workdps(40):
                    for i in np.flatnonzero(~np.isfinite(tail)):
                        shape=mp.mpf(float(s[i]));yy=mp.mpf(y)
                        tail[i]=float(2*yy**(shape/2)*mp.besselk(shape,2*mp.sqrt(yy))/mp.gamma(shape))
        else:tail=gammaincc(s,y)
        return float(tail[0]),float(weights@(1-tail))

    def optimum(self,alpha=.01):
        if not 0<alpha<1:raise ValueError('Interior blank limit required.')
        y=brentq(lambda y:self.errors(y)[0]-alpha,1e-10,100.)
        return dict(y=y,blank=alpha,miss=self.errors(y)[1],evidence='numerical limiting law; not finite-h inference')