"""Paper-conditional mission arithmetic, with fresh outward interval evaluation.

The process probability theorem is imported, not re-proved by this code.
"""
from fractions import Fraction as Q
from math import ceil,factorial,comb
from mpmath import mp,iv
import sympy as sp
from chemistry import ceildiv,restart,dot,I
mp.dps=70;iv.dps=70
KAPPA=Q(1839,8750000000000)


def num(ctx,x):
    x=Q(str(x));return ctx.mpf(x.numerator)/x.denominator


def exp_lower(x,terms=160):
    x=Q(x)
    if x<0:raise ValueError('Positive-exponential Taylor bound requires x >= 0.')
    term=total=Q(1)
    for k in range(1,terms):term*=x/k;total+=term
    return total


class MissionCertificate:
    def __init__(self,V,reactor,bath,refined=True,sharp=False):
        if type(V) is not int or V<10**6 or not reactor.admitted(bath):raise ValueError('One-cycle theorem requires integer V >= 1e6 and the admitted finite-bath class.')
        if sharp and (not refined or bath.inventory!=bath.capacity or reactor.d!=Q(1,50)):raise ValueError('Halved service allowance requires refined pure-total bath M=R and d=1/50.')
        self.V=V;self.refined=refined;self.sharp=sharp;self.rate=KAPPA if refined else Q(1,10**10)

    def logs(self,ctx=mp):
        V=ctx.mpf(self.V);ln=ctx.log
        terms={'pulse_stock':-num(ctx,'1177/1000000000')*V,'pulse_material':ln(4)-V/100000,
            'free_counter':-V/40000000,'phase_occupation':ln(100)-num(ctx,self.rate)*V,'free_clock_early':-V,'free_clock_collection':-V/200,
            'template_counter':-V/100000,'foods':ln(2)-V/300,'gross_service':-V/2000,
            'material_initial_twice':ln(8)-V/2000,'material_leak_twice':ln(96000)+ln(V)-V/2000+num(ctx,'1/50'),
            'recovery':-3*V/5000000,'recovery_clock':-V/2,'residence_initial':-V/10000,'residence_leak':ln(12000)+ln(V)-V/10000+num(ctx,'9/500'),
            'terminal_material':ln(4)-V/320000,'terminal_clock':ln(4)-V}
        if self.sharp:terms['additional_sharp_service']=-V/2000
        return terms

    def log_error(self,ctx=mp):
        logs=self.logs(ctx);shift=logs['phase_occupation'];return shift+ctx.log(sum(ctx.exp(x-shift) for x in logs.values()))

    def evaluate(self,m):
        if type(m) is not int or m<1:raise ValueError('Positive integer mission length required.')
        L=self.log_error(iv);e=iv.exp(L);failure=m*e
        lower=(1-e)**m if bool(e.b<iv.mpf(1).a) else iv.mpf(0)
        digits=max(0,int(mp.floor(mp.mpf(lower.a)*10**15)))
        while digits and not bool(num(iv,Q(digits,10**15)).b<=lower.a):digits-=1
        return dict(V=self.V,cycles=m,refined_phase=self.refined,sharp_service=self.sharp,
            single_cycle_log_interval=str(L),mission_failure_upper_display=mp.nstr(mp.mpf(failure.b),28),
            joint_product_lower_rational=str(Q(digits,10**15)),joint_product_lower_display=str(Q(digits,10**15).numerator/Q(digits,10**15).denominator),
            component_log10={k:mp.nstr(v/mp.log(10),22) for k,v in self.logs().items()},
            scope='Initial restart, admitted pulse and fixed-rate source required. Conditional iteration on actual returns, not independent cycles. Full error sum remains untruncated; probability lower bound is clipped at zero only if vacuous.')

    def inventory(self,m,initial_I=None):
        if m<1:raise ValueError('At least one cycle required.')
        V=self.V;bound=161*V//160 if initial_I is None else initial_I
        return dict(collected_I=m*ceildiv(V,56),included_collected_X=m*ceildiv(V,1080),each_food_allowance=5*m*V,gross_service_allowance=m*(V//(10 if self.sharp else 5)),final_I_floor=ceildiv(V,28),initial_I_upper=bound,net_synthesis_lower=m*ceildiv(V,56)+ceildiv(V,28)-bound)


def design(m,failure,rho,template_demand=0,free_demand=0):
    failure=Q(failure);rho=Q(rho)
    if type(m) is not int or m<1 or not 0<failure<1 or not 0<rho<1 or any(type(v) is not int or v<0 for v in (template_demand,free_demand)):raise ValueError('Invalid mission design.')
    term=iv.log(101*m/num(iv,failure))/num(iv,KAPPA)
    V=max(5*10**10,int(mp.ceil(mp.mpf(term.b))),ceildiv(56*template_demand,m),ceildiv(1080*free_demand,m))
    target=Q(101*m)/failure;value=exp_lower(KAPPA*V)
    if value<target:raise ArithmeticError('Taylor degree insufficient to certify the proposed scale.')
    gross=m*(V//10);R=max(1,ceil(Q(gross)/rho))
    return dict(V=V,R=R,cycles=m,failure=str(failure),rho=str(rho),requested_I=template_demand,requested_X=free_demand,taylor_lower=str(value),taylor_target=str(target),gross_service_allowance=gross,
        food_operation_each=5*m*V,food_with_one_spare_each=5*m*V+1,scope='Sufficient pure-bath design at d=1/50, r=20, established all-free-X restart. Preparation food charge must be added separately.')


def loaded_force_capacity(gross,zeta):
    """Outward-rounded sufficient R; caller supplies the LOADED gross allowance."""
    if gross<0 or Q(zeta)<=0:raise ValueError('Nonnegative gross bound and positive force tolerance required.')
    z=num(iv,zeta);e=iv.exp(z);bound=iv.mpf(gross)*(e+1)/(e-1)
    return max(1,int(mp.ceil(mp.mpf(bound.b))))


def bath_energy(R,J,loaded=False):
    if type(R) is not int or R<1 or type(J) is not int or not (-R<=J<=R if loaded else 0<=J<=R):raise ValueError('Endpoint current outside the bath support.')
    # Symbolic factorials remain unevaluated, including at enormous copy scale.
    if loaded:mix=sp.log(sp.factorial(R-J,evaluate=False))+sp.log(sp.factorial(R+J,evaluate=False))-2*sp.log(sp.factorial(R,evaluate=False))
    else:mix=-sp.log(sp.binomial(R,J,evaluate=False))
    exact=-J*sp.log(80000000000)+mix
    # High precision log-gamma is only a displayed approximation, never the identity.
    mixing=(mp.loggamma(R-J+1)+mp.loggamma(R+J+1)-2*mp.loggamma(R+1)) if loaded else (mp.loggamma(R-J+1)+mp.loggamma(J+1)-mp.loggamma(R+1))
    value=-J*mp.log(80000000000)+mixing
    return dict(R=R,J=J,loaded=loaded,delta_G_over_kBT=str(exact),delta_G_display=mp.nstr(value,28),bath_free_energy_decrease_display=mp.nstr(-value,28),scope='Endpoint state-function change, not gross turnover or extracted useful work; apparatus costs excluded.')


def neighbour_balance(N,V,bath,reactor):
    if reactor.d==0 or N[2]<1 or bath.fuel<1:raise ValueError('Supported active forward step required.')
    dest=tuple(n+j for n,j in zip(N,reactor.internal.channels[18].jump));after=bath.step(18)
    ratio=reactor.propensities(N,V,bath)[18]/reactor.propensities(dest,V,after)[19]
    expected=Q(8000000000*bath.fuel*V*N[2],(bath.waste+1)*(N[0]+1)*(N[1]+1))
    if ratio!=expected:raise ArithmeticError('Neighbour-state detailed balance failed.')
    return dict(propensity_ratio=str(ratio),exp_negative_free_energy_change=str(expected),post_forward_waste=after.waste)