"""Exact bounded-integer lower semigroup for the pure productive core.

Adapted from the manuscript check_paper.py; see provenance.json.
"""
from fractions import Fraction as F
from math import comb,factorial,ceil
import numpy as np
from scipy.sparse import coo_matrix


def binomial_event(n,lo,hi,p):
    p=F(p)
    if not isinstance(n,int) or n<0 or not 0<=p<=1:raise ValueError('Invalid binomial source.')
    lo=max(0,lo);hi=min(n,hi)
    if hi<lo:return F(0)
    a,d=p.numerator,p.denominator
    return F(sum(comb(n,j)*a**j*(d-a)**(n-j) for j in range(lo,hi+1)),d**n)


class PureCore:
    def __init__(self,inventory=80,minimum=8):
        if not isinstance(inventory,int) or not isinstance(minimum,int) or not 1<=minimum<=inventory:raise ValueError('Positive integer core/threshold required.')
        self.K=inventory;self.minimum=minimum
        self.states=tuple((n,p) for n in range(1,inventory+1) for p in range(inventory-n+1));self.index={s:i for i,s in enumerate(self.states)}
        entries=[];exits=[]
        for i,(n,p) in enumerate(self.states):
            f=inventory-n-p;total=0
            for dst,rate in [((n+1,p),100*n*f),((n-1,p),n*(n-1)),((n,p+1),100*n*f),((n,p-1),n*p)]:
                if rate:
                    if dst not in self.index:raise ArithmeticError('Core transition exits the exact state space.')
                    entries.append((i,self.index[dst],rate));total+=rate
            exits.append(total)
        self.lam=ceil(F(max(exits),100))
        if self.lam==0:raise ValueError('This uniformization implementation requires a nonzero generator.')
        self.D=100*self.lam;entries.extend((i,i,self.D-rate) for i,rate in enumerate(exits))
        rr,cc,aa=zip(*entries);self.A=coo_matrix((np.array(aa,dtype=np.int64),(rr,cc)),shape=(len(self.states),)*2).tocsr()
        if np.any(self.A.data<0) or not np.all(np.asarray(self.A.sum(axis=1)).ravel()==self.D):raise ArithmeticError('Uniformization rows are not stochastic.')
    def payoff(self,partition=F(1,2),recovery=F(1)):
        split={n:binomial_event(n,self.minimum,n-self.minimum,partition) for n in range(1,self.K+1)}
        collect={p:binomial_event(p,4,p,recovery) for p in range(self.K+1)}
        return [split[n]*collect[p] for n,p in self.states]
    def certify(self,partition=F(1,2),recovery=F(1),duration=F(1),scale=2**31):
        steps=self.lam*F(duration)
        if steps.denominator!=1 or steps<=0 or not isinstance(scale,int) or scale<1:raise ValueError('Positive integer number of 1/lambda blocks and positive scale required.')
        e_minus=sum((F((-1)**j,factorial(j)) for j in range(26)),F(0));weights=[int(scale*e_minus/factorial(j)) for j in range(19)]
        if self.D*scale>=2**63 or scale*sum(weights)>=2**63 or sum(weights)>scale:raise OverflowError('Requested grid violates signed-64-bit accumulator proof.')
        values=np.array([int(scale*p) for p in self.payoff(partition,recovery)],dtype=np.int64);last=max(j for j,w in enumerate(weights) if w)
        for _ in range(int(steps)):
            term=values;acc=weights[0]*term
            for j in range(1,last+1):term=self.A.dot(term)//self.D;acc+=weights[j]*term
            values=acc//scale
        if values.min()<0 or values.max()>scale:raise ArithmeticError('Lower iterate left its invariant interval.')
        starts=[self.index[(n,0)] for n in range(self.minimum,self.K+1)];worst=min(starts,key=lambda i:int(values[i]))
        return dict(lower=F(int(values[worst]),scale),minimum_state=self.states[worst],values=values,scale=scale,lambda_=self.lam,denominator=self.D,
            state_count=len(self.states),blocks=int(steps),last_nonzero_weight=last,matrix_accumulator_bound=self.D*scale,weighted_accumulator_bound=scale*sum(weights),weight_sum=sum(weights),
            scope='Exact downward integer recurrence, not a floating-point exponential. Stochastic-semigroup soundness uses the manuscript theorem; no Lean execution is claimed.')


def cluster_partition(weights,minimum=8,probability=F(1,2)):
    p=F(probability)
    if not 0<=p<=1 or any(not isinstance(w,int) or w<1 for w in weights):raise ValueError('Positive integer clusters and a valid partition probability required.')
    coefficients=[F(1)]
    for weight in weights:
        updated=[F(0)]*(len(coefficients)+weight)
        for j,a in enumerate(coefficients):updated[j]+=a*(1-p);updated[j+weight]+=a*p
        coefficients=updated
    if sum(coefficients)!=1:raise ArithmeticError('Partition polynomial lost probability mass.')
    return sum(coefficients[minimum:sum(weights)-minimum+1],F(0)) if sum(weights)>=2*minimum else F(0)