#!/usr/bin/env python3
"""Exact rational budgets and 80-place outward fixed-point evaluation.
No external dependencies. Finite checks supplement the general proofs in the pages.
"""
import argparse
from dataclasses import dataclass
from fractions import Fraction as Q
from functools import lru_cache
from math import factorial,isqrt
from pathlib import Path
import json

DIGITS=80
SCALE=10**DIGITS
CHECKS=0

def check(ok,label):
    global CHECKS
    CHECKS+=1
    if not ok:raise ArithmeticError(label)

def floorq(q):return q.numerator//q.denominator

def ceilq(q):return -((-q.numerator)//q.denominator)

@dataclass(frozen=True)
class Interval:
    lo:int
    hi:int
    def __post_init__(self):
        if self.lo>self.hi:raise ValueError('reversed interval')
    @staticmethod
    def rational(x):
        q=Q(x)*SCALE
        return Interval(floorq(q),ceilq(q))
    def __add__(self,other):
        other=as_interval(other)
        return Interval(self.lo+other.lo,self.hi+other.hi)
    __radd__=__add__
    def __neg__(self):return Interval(-self.hi,-self.lo)
    def __sub__(self,other):return self+-as_interval(other)
    def __rsub__(self,other):return as_interval(other)+-self
    def __mul__(self,other):
        other=as_interval(other)
        p=[a*b for a in [self.lo,self.hi] for b in [other.lo,other.hi]]
        return Interval(min(p)//SCALE,-((-max(p))//SCALE))
    __rmul__=__mul__
    def reciprocal(self):
        if self.lo<=0<=self.hi:raise ZeroDivisionError('interval contains zero')
        return Interval(floorq(Q(SCALE*SCALE,self.hi)),ceilq(Q(SCALE*SCALE,self.lo)))
    def __truediv__(self,other):return self*as_interval(other).reciprocal()
    def __rtruediv__(self,other):return as_interval(other)*self.reciprocal()
    def __pow__(self,n):
        if not isinstance(n,int) or n<0:raise ValueError('nonnegative integer exponent')
        result=Interval.rational(1)
        for _ in range(n):result=result*self
        return result
    def abs_upper(self):return max(abs(self.lo),abs(self.hi))
    def absolute(self):
        if self.lo>=0:return self
        if self.hi<=0:return -self
        return Interval(0,self.abs_upper())
    def contains(self,q):return self.lo<=Q(q)*SCALE<=self.hi
    def intersects(self,other):return self.lo<=other.hi and other.lo<=self.hi
    def record(self):return {'lower':decimal(self.lo),'upper':decimal(self.hi)}

def as_interval(x):return x if isinstance(x,Interval) else Interval.rational(x)
def decimal(n):
    sign='-' if n<0 else '';a,b=divmod(abs(n),SCALE)
    return sign+str(a)+'.'+str(b).zfill(DIGITS)

def atan_reciprocal(m,terms):
    # Alternating series, decreasing magnitudes; adjacent sums bracket atan(1/m).
    a=sum((Q(1,(2*j+1)*m**(2*j+1)) if j%2==0 else -Q(1,(2*j+1)*m**(2*j+1))) for j in range(terms))
    next_term=Q((-1)**terms,(2*terms+1)*m**(2*terms+1))
    return Interval(floorq(min(a,a+next_term)*SCALE),ceilq(max(a,a+next_term)*SCALE))

PI=16*atan_reciprocal(5,100)-4*atan_reciprocal(239,24)
ROOT2=Interval(isqrt(2*SCALE*SCALE),isqrt(2*SCALE*SCALE)+1)

@lru_cache(maxsize=None)
def sinpi_reduced(r):
    if not -Q(1,2)<=r<=Q(1,2):raise ValueError('reduced argument')
    if r==0:return Interval.rational(0)
    x=PI*r
    check(x.abs_upper()<2*SCALE,'Taylor reduced argument')
    x2=x*x;term=x;answer=x
    for j in range(1,40):
        term=-term*x2/Q((2*j)*(2*j+1));answer=answer+term
    # Polynomial degree 79. Degree 80 has the same value; Lagrange remainder
    # <= |x|^81/81! <= 2^81/81!, independently of alternating signs.
    e=ceilq(Q(2**81,factorial(81))*SCALE)
    return Interval(answer.lo-e,answer.hi+e)

def sinpi(x):
    x=Q(x);k=floorq(x+Q(1,2));r=x-k
    ans=sinpi_reduced(r)
    return -ans if k%2 else ans

@lru_cache(maxsize=None)
def sinc(x):
    x=Q(x)
    if not x:return Interval.rational(1)
    if x.denominator==1:return Interval.rational(0)
    return sinpi(x)/(PI*x)

def phi(x,delta=Q(1,2)):return sinc(Q(x))*sinc(delta*Q(x))

def signal(t):return sinc(Q(t)/4)**2

def sample(k):
    if k==0:return Interval.rational(1)
    if k%4==0:return Interval.rational(0)
    return (8 if k%2 else 16)/(PI**2*k*k)

def critical_sample(k):
    if k==0:return Interval.rational(1)
    if k%2==0:return Interval.rational(0)
    return 4/(PI**2*k*k)

def critical_tail_squared(N):
    if type(N) is not int or N<1:raise ValueError('positive truncation index')
    return 32/(3*PI**4*N**3)

def localized_tail(N):
    if type(N) is not int or N<4:raise ValueError('N>=4 for the declared observation window')
    return 256/(3*PI**4*N**3)

def general_budget(B,h,T,N,M,eta):
    B,h,T,M,eta=map(Q,(B,h,T,M,eta))
    if B<=0 or h<=0 or h*B>=Q(1,2) or T<0 or type(N) is not int or N<=T/h or M<0 or eta<0:
        raise ValueError('invalid bandwidth, margin, window, bound, or noise contract')
    delta=1-2*h*B
    tail=2*M/(PI**2*delta*(N-T/h))
    noise=eta*(1+12/(PI**2*delta))
    return tail,noise

def critical_partial(t,N):return sum((critical_sample(k)*sinc(Q(t)/2-k) for k in range(-N,N+1)),Interval.rational(0))
def localized_partial(t,N):return sum((sample(k)*phi(Q(t)-k) for k in range(-N,N+1)),Interval.rational(0))

def half_point_exact_partial(N):
    # phi(1/2-k)=sigma_k sqrt(2)/(pi²(1/2-k)²), sigma=+,+,-,- modulo 4.
    A=Q(0);B=Q(0)
    for k in range(-N,N+1):
        sign=1 if k%4 in [0,1] else -1;r=Q(1,2)-k
        if k==0:A+=sign/r**2
        elif k%4:B+=Q(sign*(8 if k%2 else 16),k*k)/r**2
    return ROOT2*(A/PI**2+B/PI**4),A,B

def main():
    global CHECKS
    CHECKS = 0
    sinpi_reduced.cache_clear()
    sinc.cache_clear()
    check(PI.lo>3*SCALE and PI.hi<4*SCALE,'Machin enclosure coarse range')
    check(ROOT2.lo**2<=2*SCALE*SCALE<ROOT2.hi**2,'integer square root enclosure')
    # Exact integral of [4(1-4|xi|)]² on [-1/4,1/4].
    energy=32*(Q(1,4)-4*Q(1,4)**2+Q(16,3)*Q(1,4)**3)
    check(energy==Q(8,3),'triangular spectrum energy')
    for k in range(-120,121):
        check(signal(k).intersects(sample(k)),'all three integer sample classes')
        check(signal(2*k).intersects(critical_sample(k)),'critical odd/even samples')
        check(phi(k).contains(1 if k==0 else 0),'cardinal kernel')
    tolerance=Q(1,1000)
    critical_first=None
    for N in range(1,49):
        bound=critical_tail_squared(N)
        if bound.hi<=tolerance*tolerance*SCALE:critical_first=N;break
        check(bound.lo>tolerance*tolerance*SCALE,'failed critical bound certified')
    check(critical_first==48,'first certified critical N')
    exact_one=8/PI**2
    for N in [1,2,4,8,16,32,48,64]:
        partial=critical_partial(1,N);err=(exact_one-partial).abs_upper()
        bound=critical_tail_squared(N)
        check(err*err<=bound.lo*SCALE,'critical point error versus global theorem bound')
        normpartial=2*sum((critical_sample(k)**2 for k in range(-N,N+1)),Interval.rational(0))
        remainder=as_interval(energy)-normpartial
        check(remainder.lo>=0 and remainder.hi<=2*bound.lo,'critical Parseval and L2 tail scale')
    # True L2 null-sample signal at h=3: two shifted intervals of half-width 1/20.
    low=Q(1,6)-Q(1,20);high=Q(1,6)+Q(1,20)
    check((low,high)==(Q(7,60),Q(13,60)) and high<Q(1,4),'undersampling spectral endpoints')
    for k in range(-100,101):check((sinpi(k)*sinc(Q(3*k,10))).contains(0),'all tested coarse samples vanish')
    distinction=sinc(Q(3,20));check(distinction.lo>0,'nonzero indistinguishable input at 3/2')
    # Exact trapezoid geometry, including the edge touching adjacent periodic support.
    for b in [Q(1,10),Q(1,4),Q(2,5),Q(49,100)]:
        delta=1-2*b
        for j in range(-200,201):
            w=Q(j,100)
            overlap=max(Q(0),min(Q(1,2),w+delta/2)-max(-Q(1,2),w-delta/2))
            W=overlap/delta
            expected=Q(1) if abs(w)<=b else (1-b-abs(w))/delta if abs(w)<1-b else Q(0)
            check(W==expected,'trapezoid from exact overlap')
            if abs(w)<=b:check(W==1,'exact passband')
            for shift in [-2,-1,1,2]:
                if abs(w-shift)<=b:check(W==0,'adjacent copies killed including endpoints')
    eta=Q(1,100000);noise=eta*(1+24/PI**2);tol=Q(1,10000)
    first=None
    for N in range(4,25):
        total=localized_tail(N)+noise
        if total.hi<=tol*SCALE:first=N;break
        check(total.lo>tol*SCALE,'failed noisy budget certified')
    check(first==24,'first certified noisy window N')
    for N in [4,8,12,24,32]:
        bound=localized_tail(N)
        for j in range(-16,17):
            t=Q(j,8);err=(signal(t)-localized_partial(t,N)).abs_upper()
            check(err<=bound.lo,'sampled rational observation points inside analytic bound')
    partial,A,B=half_point_exact_partial(24);generic=localized_partial(Q(1,2),24)
    check(partial.intersects(generic),'independent radical and Taylor evaluation of finite sum')
    evaluation_budget=Q(1,10**70)
    check(generic.hi-generic.lo<=evaluation_budget*SCALE,'finite evaluation width budget')
    check((localized_tail(24)+noise+evaluation_budget).hi<=tol*SCALE,'tail noise and evaluation all fit')
    target=16*(2-ROOT2)/PI**2
    check(target.intersects(signal(Q(1,2))),'exact half point signal')
    check((target-partial).abs_upper()<=localized_tail(24).lo,'half point exact sum and tail')
    rows=[];previous=Interval.rational(0)
    for N in [4,8,16,32,64,128,256]:
        harmonic=sum((1/abs(Q(k)-Q(1,2)) for k in range(-N,N+1)),Q(0))
        raw=eta*harmonic/PI
        window_sum=sum((phi(Q(1,2)-k).absolute() for k in range(-N,N+1)),Interval.rational(0))
        # No zero crossing occurs for the half-integer kernels; compare upper to lower.
        check(raw.lo>previous.hi,'adversarial sinc amplification strictly increases')
        check((eta*window_sum).hi<=noise.lo,'bounded-noise localized certificate')
        previous=raw;rows.append({'N':N,'sinc_noise':raw.record(),'localized_adversarial_noise':(eta*window_sum).record()})
    rejects=0
    invalid=[(Q(3,4),1,2,24,1,eta),(Q(1,2),1,2,24,1,eta),(Q(1,4),0,2,24,1,eta),(0,1,2,24,1,eta),(Q(1,4),1,-1,24,1,eta),(Q(1,4),1,2,2,1,eta),(Q(1,4),1,2,24,-1,eta),(Q(1,4),1,2,24,1,-eta),(Q(1,4),1,2,True,1,eta)]
    for args in invalid:
        try:general_budget(*args)
        except ValueError:rejects+=1
        else:raise ArithmeticError('invalid certificate accepted')
    check(rejects==len(invalid),'bandwidth/window/noise contract rejection')
    general,gn=general_budget(Q(1,4),1,2,24,1,eta)
    check(gn==noise and general.lo>localized_tail(24).hi,'generic M bound distinct from sample decay')
    return {'status':'PASS','checks':CHECKS,'arithmetic':'80-place integer outward intervals; rational Machin bounds; integer sqrt(2); reduced Taylor remainder','energy':str(energy),'critical':{'first_N':48,'N47_squared_bound':critical_tail_squared(47).record(),'N48_squared_bound':critical_tail_squared(48).record(),'target_t1':exact_one.record()},'undersampling':{'spectral_positive_endpoints':[str(low),str(high)],'g_at_3_over_2':distinction.record()},'localized':{'first_N':24,'N23_total_bound':(localized_tail(23)+noise).record(),'N24_total_bound':(localized_tail(24)+noise).record(),'noise_bound':noise.record(),'generic_tail_N24':general.record(),'sample_decay_tail_N24':localized_tail(24).record(),'half_point_partial':partial.record(),'half_point_target':target.record(),'half_point_coefficients_A_B':[str(A),str(B)],'half_point_evaluation_width_budget':str(evaluation_budget)},'adversarial_noise':rows,'invalid_contracts_rejected':rejects,'scope':'Theorems and uniform tail/noise bounds are proved in the pages. Rational point evaluations are finite diagnostics, not a substitute for whole-window proofs.'}

if __name__=='__main__':
    parser=argparse.ArgumentParser();parser.add_argument('--output',type=Path,required=True);args=parser.parse_args()
    result=main();args.output.parent.mkdir(parents=True,exist_ok=True);args.output.write_text(json.dumps(result,ensure_ascii=False,indent=2)+'\n');print(json.dumps({'status':result['status'],'checks':result['checks'],'output':str(args.output)}))