#!/usr/bin/env python3
"""Exact rational and outward-interval checks for equivalence/content contracts.
Only Python standard library. --output chooses the result path.
No black-box CDF, random simulation, or floating-point acceptance decisions.
"""
from fractions import Fraction as Q
from math import isqrt,comb,factorial
from collections import Counter
from itertools import product
from pathlib import Path
import argparse,json
BITS=220;SCALE=1<<BITS;checks=Counter()
def check(ok,label):
 checks[label]+=1
 if not ok:raise ArithmeticError(label)
def down(q):return Q(q.numerator*SCALE//q.denominator,SCALE)
def up(q):return -down(-q)
class I:
 def __init__(self,lo,hi=None):
  lo=Q(lo);hi=lo if hi is None else Q(hi)
  if lo>hi:raise ValueError('reversed interval')
  self.lo=down(lo);self.hi=up(hi)
 def __add__(self,b):
  b=asI(b);return I(self.lo+b.lo,self.hi+b.hi)
 __radd__=__add__
 def __neg__(self):return I(-self.hi,-self.lo)
 def __sub__(self,b):return self+-asI(b)
 def __rsub__(self,b):return asI(b)+-self
 def __mul__(self,b):
  b=asI(b);v=[x*y for x,y in product([self.lo,self.hi],[b.lo,b.hi])];return I(min(v),max(v))
 __rmul__=__mul__
 def inv(self):
  if self.lo<=0<=self.hi:raise ValueError('division interval meets zero')
  return I(1/self.hi,1/self.lo)
 def __truediv__(self,b):return self*asI(b).inv()
 def __rtruediv__(self,b):return asI(b)*self.inv()
 def square(self):
  return I(0 if self.lo<=0<=self.hi else min(self.lo*self.lo,self.hi*self.hi),max(self.lo*self.lo,self.hi*self.hi))
 def sqrt(self):
  if self.lo<0:raise ValueError('negative square root')
  a=isqrt(self.lo.numerator*SCALE*SCALE//self.lo.denominator)
  b=isqrt(self.hi.numerator*SCALE*SCALE//self.hi.denominator)
  if Q(b*b,SCALE*SCALE)<self.hi:b+=1
  return I(Q(a,SCALE),Q(b,SCALE))
 def contains(self,q):return self.lo<=q<=self.hi
 def width(self):return self.hi-self.lo
 def show(self,digits=18):
  base=10**digits
  def fixed(n):
   sign='-' if n<0 else '';n=abs(n);return sign+str(n//base)+'.'+str(n%base).zfill(digits)
  return {'lo':fixed(self.lo.numerator*base//self.lo.denominator),'hi':fixed(-((-self.hi.numerator*base)//self.hi.denominator))}
def asI(x):return x if isinstance(x,I) else I(x)
def atan_recip(q,N):
 a=Q(0)
 for j in range(N):a+=Q((-1)**j,(2*j+1)*q**(2*j+1))
 nxt=Q((-1)**N,(2*N+1)*q**(2*N+1));return I(min(a,a+nxt),max(a,a+nxt))
PI=16*atan_recip(5,100)-4*atan_recip(239,35)
ROOT2PI=(2*PI).sqrt()
def exp_point(x):
 x=Q(x)
 if x<0:return exp_point(-x).inv()
 z=x;h=0
 while z>Q(1,8):z/=2;h+=1
 term=Q(1);acc=term;N=0
 while True:
  nxt=term*z/(N+1);ratio=z/(N+2)
  rem=nxt/(1-ratio)
  if rem<Q(1,1<<200):break
  N+=1;term=nxt;acc+=term
 ans=I(acc,acc+rem)
 for _ in range(h):ans=ans.square()
 return ans
def expI(a):
 a=asI(a);return I(exp_point(a.lo).lo,exp_point(a.hi).hi)
def phi(a):a=asI(a);return expI(-a.square()/2)/ROOT2PI
def Phi_point(x):
 x=Q(x)
 if x<0:return 1-Phi_point(-x)
 term=x;acc=term;j=0
 while True:
  nxt=term*x*x*(2*j+1)/(2*(j+1)*(2*j+3))
  # From this point onward term ratios are <1, hence alternating remainder.
  if j+1>x*x/2 and nxt<Q(1,1<<195):break
  j+=1;term=nxt;acc+=(-1)**j*term
 signed=(-1)**(j+1)*nxt
 return I(Q(1,2))+I(min(acc,acc+signed),max(acc,acc+signed))/ROOT2PI
def Phi(a):
 a=asI(a);return I(Phi_point(a.lo).lo,Phi_point(a.hi).hi)
def nct4(t,delta):
 t=asI(t);delta=asI(delta);A=4+t.square();root=A.sqrt();b=t*delta/A;q=t*delta/root;K=expI(-2*delta.square()/A)/root
 return Phi(-delta)+t*K*((1+2*b.square()+2/A)*Phi(q)+2*b/root*phi(q))
# Elementary bound checks are independent exact identities, not floating-point tests.
for q in [Q(j,17) for j in range(200)]:
 z=I(q).sqrt();check(z.lo*z.lo<=q<=z.hi*z.hi,'sqrt_outward_square')
for j in range(-20,21):
 q=Q(j,7);v=exp_point(q)*exp_point(-q);check(v.contains(1),'exp_reciprocal_identity')
 check((Phi_point(q)+Phi_point(-q)).contains(1),'normal_reflection_identity')
check(PI.lo>Q(3141592653589793238,10**18) and PI.hi<Q(3141592653589793239,10**18),'pi_decimal_bracket')
# Exact finite distribution: ten Bernoulli trials, two outside-null components.
def binmass(n,p,j):return comb(n,j)*p**j*(1-p)**(n-j)
def tail_ge(n,p,j):return sum(binmass(n,p,k) for k in range(j,n+1))
def tail_le(n,p,j):return sum(binmass(n,p,k) for k in range(j+1))
accepted=[]
for x in range(11):
 pL=tail_ge(10,Q(1,5),x);pU=tail_le(10,Q(4,5),x)
 check(pL==tail_le(10,Q(4,5),10-x),'binomial_reflection')
 if max(pL,pU)<Q(1,20):accepted.append(x)
check(accepted==[5],'binomial_exact_rejection_set')
check(tail_ge(10,Q(1,5),5)==Q(320249,9765625),'binomial_component_p')
size=binmass(10,Q(1,5),5);check(size==Q(258048,9765625)<Q(1,20),'binomial_size')
for j in range(101):
 p=Q(j,100)
 if p<=Q(1,5) or p>=Q(4,5):check(binmass(10,p,5)<=size,'binomial_size_grid_crosscheck')
# Coupled rejection indicators with the same marginals: 2x2 table enumeration.
for p,q in product([Q(j,20) for j in range(21)],repeat=2):
 lo=max(Q(0),p+q-1);hi=min(p,q)
 # Deduplicate degenerate Frechet segments at marginal probabilities zero/one.
 for joint in sorted({lo+Q(j,20)*(hi-lo) for j in range(21)}):
  cells=[joint,p-joint,q-joint,1-p-q+joint]
  check(min(cells)>=0 and sum(cells)==1,'valid_dependent_rejection_table')
  check(joint<=min(p,q),'IUT_intersection_bound')
  if p<=Q(1,20) or q<=Q(1,20):check(joint<=Q(1,20),'union_null_global_level')
check(1-(1-Q(1,20))**2==Q(39,400)>Q(1,20),'OR_independent_error_inflation')
# Projection nuisance direction never enters the numerator or denominator.
m=[Q(1),Q(1),Q(10),Q(0),Q(0),Q(0),Q(0)]
check(sum(v*v for v in m[:2])==2 and sum(v*v for v in m)==102,'projected_noncentrality')
# Independent Poisson/Beta finite sum and closed expression, all outward.
x=Q(3);b=x/(x+2);M=24;w=exp_point(-1);partial=I(0)
for j in range(M+1):
 beta=b**(j+1)*(j+2-(j+1)*b)
 check(0<=beta<=1,'beta_component_cdf_range')
 partial+=w*beta;w=w/(j+1)  # next weight after j
# At j=0 this is w1=w0; at j=M it is w_{M+1}.
tail=w/(1-Q(1,M+2));mixture=I(partial.lo,(partial+tail).hi)
closed=b*(2-b*b)*exp_point(b-1)
check(mixture.lo<=closed.lo<=closed.hi<=mixture.hi,'F_closed_inside_independent_mixture')
check(tail.hi<Q(1,10**24),'Poisson_omitted_mass_budget')
power=1-closed
check(power.lo>Q(3404050747009309,10**16) and power.hi<Q(3404050747009310,10**16),'F_power_decimal_bracket')
# Signed noncentral t; central identity uses an independently derived formula.
for t in [Q(j,3) for j in range(-12,13)]:
 T=I(t);h=T/(4+T.square()).sqrt();central=I(Q(1,2))+Q(3,4)*(h-h*h*h/3)
 v=nct4(T,0);check(max(v.lo,central.lo)<=min(v.hi,central.hi),'central_t4_integrated_density_identity')
 for delta in [Q(-2),Q(0),Q(1),Q(2)]:
  v=nct4(t,delta);mirror=nct4(-t,-delta)
  check((v+mirror).contains(1),'signed_nct_reflection')
  check(v.lo>=0 and v.hi<=1,'nct_cdf_range')
root5=I(5).sqrt();ptost=(27-11*root5)/54
check(ptost.hi<Q(1,20),'TOST_n5_certified')
check((1-nct4(root5,0)-ptost).contains(0),'TOST_exact_formula')
# A tolerance factor bracket is a pair of strict CDF comparisons, not a root guess.
klo=Q(2809,1000);khi=Q(281,100)
Flo=nct4(root5*klo,root5);Fhi=nct4(root5*khi,root5)
check(Flo.hi<Q(949994,10**6)<Q(19,20),'tolerance_lower_factor_rejected')
check(Fhi.lo>Q(950051,10**6)>Q(19,20),'tolerance_upper_factor_certified')
# Equal-tail 90% mean CI versus 90% future prediction interval on the same data.
c_lo=Q(2131846,10**6);c_hi=Q(2131847,10**6)
check(nct4(c_lo,0).hi<Q(19,20)<nct4(c_hi,0).lo,'central_t95_bracket')
c=I(c_lo,c_hi);mean_half=c/(10*root5);pred_half=c*I(Q(6,5)).sqrt()/10
check(mean_half.hi<Q(1,10)<pred_half.lo,'mean_and_prediction_different_scale')
# Known-scale TOST: z_.95 and n=34 versus n=35 are rigorously bracketed.
zlo=Q(1644853626,10**9);zhi=Q(1644853627,10**9)
check(Phi_point(zlo).hi<Q(19,20)<Phi_point(zhi).lo,'normal_z95_bracket')
z=I(zlo,zhi);powers={}
for n in [25,34,35]:powers[str(n)]=2*Phi(I(n).sqrt()/2-z)-1
check(powers['34'].hi<Q(4,5)<powers['35'].lo,'known_scale_minimum_sample_size_35')
# Valid unequal-tail 90% interval, but its containment rule can exceed 5%.
check(Phi_point(3).lo>Q(99,100),'z99_less_than_three')
check((1-Phi_point(7)).hi<Q(1,100),'upper_normal_tail_7')
# A conservative two-sided n=2, p=.5, confidence .81 construction.
check(Phi_point(Q(27,40)).lo>Q(3,4),'two_sided_z75_upper')
check(Phi_point(Q(329,200)).lo>Q(19,20),'two_sided_z95_upper')
check(Phi_point(Q(1,8)).hi<Q(11,20),'two_sided_z55_lower')
safe_factor=(I(Q(27,40))+I(Q(329,200))/I(2).sqrt())/Q(1,8)
check(safe_factor.hi<15 and Q(9,10)**2==Q(81,100),'two_sided_conservative_factor_15')
# Ternary bracketing cannot stall even when a probe equals the target root.
def shrink(a,b,root):
 u=(2*a+b)/3;v=(a+2*b)/3
 for q in [u,v]:
  if q<root:return q,b
  if q>root:return a,q
 raise ArithmeticError('distinct probes cannot both be the root')
for root in [Q(j,31) for j in range(32)]+[Q(1,3),Q(1,2),Q(2,3)]:
 a=Q(-1);b=Q(2)
 for j in range(16):
  aa,bb=shrink(a,b,root);check(aa<=root<=bb and bb-aa<=Q(2,3)*(b-a),'two_probe_quantile_bracket');a,b=aa,bb
result={'status':'PASS','checks':sum(checks.values()),'groups':dict(checks),'arithmetic':'220-bit outward dyadic rational intervals; alternating arctangent/normal integral and positive exponential Taylor remainders',
 'TOST_n5_p':ptost.show(),'binomial':{'rejection_set':accepted,'component_p':'320249/9765625','size':str(size)},'F_2_4_lambda2_tail_at3':power.show(),'Poisson_tail_upper':tail.show(35),
 'normal_tolerance_n5_p_Phi1_gamma95':{'factor_open_bracket':[str(klo),str(khi)],'cdf_at_lower':Flo.show(),'cdf_at_upper':Fhi.show(),'safe_upper_limit_for_S_point1':'281/1000'},
 'mean_90_halfwidth':mean_half.show(),'prediction_90_halfwidth':pred_half.show(),'known_scale_TOST_power':{n:v.show() for n,v in powers.items()},
 'scope':['Normal model and IUT theorems are proved in the pages, not inferred from finite checks.','General noncentral-t integrals and two-sided tolerance inversions have analytic algorithms in the pages; this executable certifies the displayed four-degree and finite examples.','The Poisson omitted mass, not a comparison of floating-point library outputs, encloses the F mixture.']}
ap=argparse.ArgumentParser();ap.add_argument('--output',type=Path);args=ap.parse_args();txt=json.dumps(result,ensure_ascii=False,indent=2)+'\n'
if args.output:args.output.parent.mkdir(parents=True,exist_ok=True);args.output.write_text(txt)
print(txt,end='')