#!/usr/bin/env python3
"""Continuous tomography: exact rational identities and 80-place outward budgets.
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 math import factorial,isqrt,comb
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)

def sqrt_interval(x):
    x=as_interval(x)
    if x.lo<0:raise ValueError('nonnegative square root input required')
    a=isqrt(x.lo*SCALE);b=isqrt(x.hi*SCALE)
    if b*b<x.hi*SCALE:b+=1
    return Interval(a,b)

def exp_negative(x):
    """exp(-x), x>=0. Positive Taylor bracket after range reduction."""
    x=as_interval(x)
    if x.lo<0:raise ValueError('nonnegative argument required')
    if x.hi==0:return as_interval(1)
    y=x;n=0
    while y.hi>SCALE//2:y=y/2;n+=1
    term=as_interval(1);s=term
    for j in range(1,81):term=term*y/j;s=s+term
    # For all 0<=y<=1/2, ratios after term 81 are <=1/(2*82).
    tail=Q(1,2)**81/factorial(81)/(1-Q(1,164))
    s=Interval(s.lo,s.hi+ceilq(tail*SCALE));ans=s.reciprocal()
    for _ in range(n):ans=ans*ans
    return ans

def gaussian_budget(K,delta,model='radial-gaussian'):
    K,delta=Q(K),Q(delta)
    if K<=0 or delta<0 or model!='radial-gaussian':raise ValueError('cutoff, noise, or model contract')
    tail=exp_negative(PI*K*K)/ROOT2
    noise=sqrt_interval(K)*delta
    return {'tail':tail,'noise':noise,'triangle':tail+noise,
            'orthogonal':sqrt_interval(tail*tail+noise*noise),
            'center':1-exp_negative(PI*K*K)}

def gaussian_parameters(A,center,theta):
    if len(A)!=2 or any(len(row)!=2 for row in A) or len(center)!=2 or len(theta)!=2:
        raise ValueError('two-dimensional inputs required')
    p,q=map(Q,A[0]);q2,r=map(Q,A[1]);cx,cy=map(Q,center);x,y=map(Q,theta)
    det=p*r-q*q
    if q!=q2 or p<=0 or det<=0 or x*x+y*y!=1:raise ValueError('symmetric positive definite matrix and unit normal required')
    d=p*x*x+2*q*x*y+r*y*y
    b=q*(x*x-y*y)+(r-p)*x*y
    c=p*y*y-2*q*x*y+r*x*x
    return dict(det=det,d=d,b=b,c=c,mean=cx*x+cy*y,
                width=det/c,frequency=(r*x*x-2*q*x*y+p*y*y)/det)

def direction(t):
    t=Q(t);return ((1-t*t)/(1+t*t),2*t/(1+t*t))

def layer_profile(layers,r):
    r=Q(r)
    if r<0:raise ValueError('radius must be nonnegative')
    for R,a in layers:
        if Q(R)<=0:raise ValueError('positive layer radii')
    # Convention only: boundary values are not identifiable from integrals.
    return sum((Q(a) for R,a in layers if r<Q(R)),Q(0))

def layer_projection(layers,s):
    s=Q(s);answer=as_interval(0)
    for R,a in layers:
        R,a=Q(R),Q(a)
        if R<=0:raise ValueError('positive layer radii')
        if abs(s)<R:answer+=2*a*sqrt_interval(R*R-s*s)
    return answer

def layer_compound_coefficient(layers,t):
    # A(AQ)/pi; exact integral of Q from t upward.
    t=Q(t)
    if t<0:raise ValueError('nonnegative squared radius')
    return sum((Q(a)*max(Q(R)**2-t,Q(0)) for R,a in layers),Q(0))

def poly_multiply(a,b):
    out={}
    for (i,j),v in a.items():
        for (k,l),w in b.items():out[i+k,j+l]=out.get((i+k,j+l),Q(0))+v*w
    return {e:v for e,v in out.items() if v}

def poly_evaluate(p,x,y):return sum((a*x**i*y**j for (i,j),a in p.items()),Q(0))

def reject(call):
    try:call()
    except (ValueError,ZeroDivisionError):check(True,'invalid certificate rejected');return
    raise ArithmeticError('invalid certificate accepted')

def main():
    global CHECKS
    CHECKS = 0
    check(3*SCALE<PI.lo<=PI.hi<4*SCALE,'Machin enclosure')
    matrices=0
    normals=[(Q(1),Q(0)),(Q(0),Q(1))]+[direction(Q(j,3)) for j in range(-5,6)]
    for p in range(1,5):
        for r in range(1,5):
            for j in range(-6,7):
                q=Q(j,3)
                if p*r<=q*q:continue
                matrices+=1
                for theta in normals:
                    z=gaussian_parameters([[p,q],[q,r]],[Q(1,2),-1],theta)
                    check(z['c']>0 and z['d']*z['c']-z['b']**2==z['det'],'rotated determinant')
                    check(z['frequency']==z['c']/z['det'],'one-dimensional spectrum matches slice')
                    check(Q(1,z['c'])/z['width']==1/z['det'],'mass squared includes amplitude')
                    for u,v in [(Q(0),Q(1)),(Q(1,3),Q(-2,5)),(Q(-3),Q(7,4))]:
                        left=z['d']*u*u+2*z['b']*u*v+z['c']*v*v
                        right=z['c']*(v+z['b']*u/z['c'])**2+z['width']*u*u
                        check(left==right,'complete square along original line')
    specific=gaussian_parameters([[2,1],[1,2]],[Q(1,2),-1],[Q(3,5),Q(4,5)])
    check(specific['c']==Q(26,25) and specific['mean']==-Q(1,2),'third projection geometry')
    check(specific['width']==Q(75,26) and specific['frequency']==Q(26,75),'third projection exponent')
    # Exact finite-angle null polynomials, omitting the nonzero Fourier factors (2pi i)^m.
    null_models=0
    for m in range(1,10):
        ns=[direction(Q(j,3)) for j in range(m)]
        poly={(0,0):Q(1)}
        for x,y in ns:poly=poly_multiply(poly,{(1,0):-y,(0,1):x})
        check(bool(poly),'nonzero product polynomial')
        for x,y in ns:
            for sigma in [Q(-3),Q(-1,2),Q(0),Q(2,3),Q(7)]:
                check(poly_evaluate(poly,sigma*x,sigma*y)==0,'all measured frequency lines annihilated')
        # A concrete nonzero frequency prevents a finite check of zeros becoming the only evidence.
        witness=next((Q(1),Q(k)) for k in range(1,m+3) if poly_evaluate(poly,Q(1),Q(k)))
        check(poly_evaluate(poly,*witness)!=0,'unmeasured nonzero witness');null_models+=1
    check(-2*Q(3,5)*Q(4,5)==-Q(24,25),'third-angle zero-distance coefficient of pi')
    # Complete budget intervals. No floating-point comparison determines success.
    budget_rows={}
    for K in [Q(1,2),Q(3,4),Q(1),Q(5,4),Q(3,2),Q(2),Q(3),Q(4)]:
        result=gaussian_budget(K,Q(1,1000))
        check(result['triangle'].hi-result['triangle'].lo<10**10,'triangle width below 1e-70')
        check(result['orthogonal'].hi<=result['triangle'].lo,'orthogonal budget improves strict triangle')
        budget_rows[str(K)]={k:v.record() for k,v in result.items()}
    one=gaussian_budget(1,Q(1,1000));five=gaussian_budget(Q(5,4),Q(1,1000))
    check(one['triangle'].lo>Q(1,100)*SCALE,'K=1 sufficient bound fails')
    check(five['triangle'].hi<Q(1,100)*SCALE,'K=5/4 sufficient bound passes')
    check(one['triangle'].lo>Q(3155685,10**8)*SCALE and one['triangle'].hi<Q(3155686,10**8)*SCALE,'published K=1 decimals')
    check(five['triangle'].lo>Q(633775,10**8)*SCALE and five['triangle'].hi<Q(633776,10**8)*SCALE,'published K=5/4 decimals')
    annuli=[]
    for K in [Q(1,3),Q(1),Q(5,4),Q(4),Q(100)]:
        for j in range(1,13):
            eps=K/2**j;a=Q(j,3)
            input_coefficient=2*eps*a*a
            output_coefficient=(K*K-(K-eps)**2)*a*a
            ratio=output_coefficient/input_coefficient
            check(ratio==K-eps/2 and ratio<K,'exact annulus squared norm ratio')
            if K==Q(5,4):annuli.append({'epsilon':str(eps),'ratio_squared':str(ratio)})
    check(Q(5,4)-Q(1,4)/2==Q(9,8),'explicit ratio greater than one')
    # Integer square-root bounds independently checked by integer squares.
    for n in range(81):
        for d in [1,3,7,19]:
            q=Q(n,d);z=sqrt_interval(q)
            check(Q(z.lo,SCALE)**2<=q<=Q(z.hi,SCALE)**2,'outward radical')
    # Abel coefficients: polynomial expansion versus a positive recurrence.
    coefficient_rows=[];C=Q(2)
    for n in range(81):
        expanded=sum((Q((-1)**j*2*comb(n,j),2*j+1) for j in range(n+1)),Q(0))
        if n:C*=Q(2*n,2*n+1)
        D=Q(comb(2*n+2,n+1),4**(n+1))
        check(C==expanded and C>0,'Abel polynomial coefficient')
        check(C*D==Q(1,n+1),'two half integrals give pi times primitive')
        if n<9:coefficient_rows.append({'degree':n,'C_n':str(C),'D_n':str(D)})
    layers=[(Q(2),Q(1)),(Q(1),Q(1))]
    checks_by_layer=0
    for model in [layers,[(Q(3,2),Q(2)),(Q(2,3),Q(-1))],[(Q(4),Q(-2)),(Q(2),Q(3)),(Q(1,2),Q(5))]]:
        breaks=sorted({Q(0)}|{R*R for R,a in model})
        for left,right in zip(breaks,breaks[1:]):
            for j in range(1,20):
                t=left+(right-left)*Q(j,20);h=min(t-left,right-t)/3
                slope=(layer_compound_coefficient(model,t+h)-layer_compound_coefficient(model,t-h))/(2*h)
                density=sum((a for R,a in model if t<R*R),Q(0))
                check(-slope==density,'full piecewise affine interval slope');checks_by_layer+=1
        for j in range(-20,21):
            s=Q(j,7)
            check(layer_projection(model,s)==layer_projection(model,-s),'even projection for signed layers')
    check(layer_profile(layers,Q(1,2))==2 and layer_profile(layers,Q(3,2))==1,'both recovered layers')
    check(layer_profile(layers,3)==0,'outside support')
    check(layer_projection(layers,0).contains(6),'central chord total')
    check(layer_projection(layers,1).intersects(2*sqrt_interval(3)),'first layer edge chord')
    check(layer_projection(layers,Q(3,2)).intersects(sqrt_interval(7)),'outer-only chord')
    check(layer_projection(layers,2).contains(0),'tangent chord zero')
    check((1-Q(1,2)**2)**2==Q(9,16),'polynomial radial query')
    # Rejected contracts include extra coordinates and false radial model assumptions.
    invalid=[lambda:gaussian_parameters([[2,1],[1,2]],[0,0,1],[1,0]),
             lambda:gaussian_parameters([[2,1,0],[1,2,0]],[0,0],[1,0]),
             lambda:gaussian_parameters([[1,2],[2,1]],[0,0],[1,0]),
             lambda:gaussian_parameters([[2,1],[0,2]],[0,0],[1,0]),
             lambda:gaussian_parameters([[2,1],[1,2]],[0,0],[1,1]),
             lambda:gaussian_budget(0,Q(1,1000)),lambda:gaussian_budget(1,-1),
             lambda:gaussian_budget(1,0,'arbitrary-finite-angle-image'),
             lambda:layer_profile(layers,-1),lambda:layer_projection([(0,1)],0),
             lambda:sqrt_interval(-1),lambda:exp_negative(-1)]
    for call in invalid:reject(call)
    return {'status':'PASS','checks':CHECKS,'arithmetic':'Fraction identities; 80-place outward integer fixed point',
            'gaussian_matrices':matrices,'unit_normals':len(normals),'finite_angle_null_models':null_models,
            'specific_projection':{k:str(v) for k,v in specific.items()},'budgets':budget_rows,
            'annulus_ratios_at_K_5_4':annuli,'abel_coefficients':coefficient_rows,'layer_interior_checks':checks_by_layer,
            'layer_values':{'1/2':'2','3/2':'1','3':'0'},'rejected_contracts':len(invalid),
            'scope':'Continuous full-angle analytic data. Finite algebraic checks supplement the general proofs, not a finite-angle scanner certificate.'}

if __name__=='__main__':
    p=argparse.ArgumentParser();p.add_argument('--output',type=Path,required=True);args=p.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)}))