#!/usr/bin/env python3
"""Exact rational identities and 70-digit cross-checks for integral/sum certificates.
Python 3.9+ and mpmath. No network, no installation, no files modified; JSON on stdout.
Analytic bounds are proved in the articles. mpmath checks are not directed interval proofs.
"""
from fractions import Fraction as F
from math import comb, factorial
import json
import mpmath as mp
mp.mp.dps=70
checks=0

def require(ok,message='check failed'):
 global checks
 checks+=1
 if not ok:raise AssertionError(message)
def mpq(x):return mp.mpf(x.numerator)/x.denominator if isinstance(x,F) else mp.mpf(x)
def close(a,b,tol='1e-50'):require(abs(a-b)<=mp.mpf(tol)*max(1,abs(b)),str((a,b)))
def bernoulli(n):
 b=[F(1)]
 for j in range(1,n+1):b.append(-sum(F(comb(j+1,k))*b[k] for k in range(j))/F(j+1))
 return b
def polynomial_b(n,b):return [F(comb(n,k))*b[n-k] for k in range(n+1)]
def poly_value(c,x):
 out=F(0)
 for v in reversed(c):out=out*x+v
 return out
def monomial_remainder(d,p,a,b,bnums):
 # Integrate B_2p(u) times the exact derivative of (k+u)^d on each unit cell.
 if d<2*p:return F(0)
 c=polynomial_b(2*p,bnums);q=d-2*p;total=F(0)
 for k in range(a,b):
  for j,v in enumerate(c):
   for ell in range(q+1):total+=v*comb(q,ell)*k**(q-ell)/F(j+ell+1)
 return -F(factorial(d),factorial(d-2*p)*factorial(2*p))*total
def tail_bounds(n):
 if not isinstance(n,int) or isinstance(n,bool) or n<1:raise ValueError('positive integer N required')
 prefix=sum((F(1,k*k) for k in range(1,n)),F(0))
 A=F(1,n)+F(1,2*n*n)+F(1,6*n**3)
 upper=prefix+A;lower=upper-F(1,16*n**5)
 return lower,upper

def main():
 result={'precision_decimal_digits':70,'exact_arithmetic':'fractions.Fraction','numeric_claim':'70-digit cross-checks; no certified directed-rounding quadrature'}
 B=bernoulli(20)
 require(B[:5]==[F(1),F(-1,2),F(1,6),F(0),F(-1,30)])
 polychecks=0
 for n in range(1,21):
  c=polynomial_b(n,B);prev=polynomial_b(n-1,B)
  require([k*c[k] for k in range(1,len(c))]==[n*x for x in prev]);polychecks+=1
  require(sum(v/F(k+1) for k,v in enumerate(c))==0);polychecks+=1
  for x in [F(j,7) for j in range(-7,15)]:
   require(poly_value(c,x+1)-poly_value(c,x)==n*x**(n-1));polychecks+=1
   require(poly_value(c,1-x)==(-1)**n*poly_value(c,x));polychecks+=1
 for j in range(1001):
  u=F(j,1000);c=polynomial_b(4,B);v=poly_value(c,u)
  require(v==u*u*(1-u)**2-F(1,30));require(F(-1,30)<=v<=F(7,240))
  require(F(-1,16)<=B[4]-v<=0)
 emcases=0
 for d in range(13):
  for p in range(1,7):
   for a,b in [(-3,-1),(-2,3),(0,1),(0,5),(1,6)]:
    total=sum(F(k)**d for k in range(a,b+1))
    approx=F(b**(d+1)-a**(d+1),d+1)+F(a**d+b**d,2)
    for j in range(1,p+1):
     order=2*j-1
     if order<=d:
      approx+=B[2*j]/factorial(2*j)*F(factorial(d),factorial(d-order))*(b**(d-order)-a**(d-order))
    rem=monomial_remainder(d,p,a,b,B)
    require(total==approx+rem);emcases+=1
 require(monomial_remainder(4,1,0,1,B)==F(-1,30))
 require(monomial_remainder(4,2,0,1,B)==0)
 # Mutation: wrong sign in the exact fourth-power remainder cannot pass.
 wrong=F(31,30)-monomial_remainder(4,1,0,1,B)
 require(wrong!=1)
 result['bernoulli_identity_checks']=polychecks
 result['finite_euler_maclaurin_exact_cases']=emcases
 result['B0_through_B20']=[str(x) for x in B]
 # 1000 independent numerical checks of exact rational tail enclosures.
 target=mp.zeta(2);prefix=F(0);boundcases=0
 for n in range(1,1001):
  A=F(1,n)+F(1,2*n*n)+F(1,6*n**3);upper=prefix+A;lower=upper-F(1,16*n**5)
  require(mpq(lower)<=target<=mpq(upper));require(upper-lower==F(1,16*n**5));boundcases+=1
  if n in [20,100,200]:
   result['zeta2_enclosure_'+str(n)]={'lower_exact':str(lower),'upper_exact':str(upper),'width_exact':str(upper-lower),'lower_decimal_display_only':mp.nstr(mpq(lower),40),'upper_decimal_display_only':mp.nstr(mpq(upper),40)}
  prefix+=F(1,n*n)
 result['tail_enclosure_N_through']=boundcases
 require(target < mpq(sum((F(1,k*k) for k in range(1,100)),F(0))+F(1,100)+F(1,2*100**2)+F(1,6*100**3)))
 for bad in [0,-1,True,1.5]:
  try:tail_bounds(bad)
  except ValueError:require(True)
  else:raise AssertionError('bad N accepted')
 # Endpoint regularization: t=u^8/2 on each half of the beta integral.
 betacases=0
 params=[mp.mpf(x) for x in ['.25','.5','.75','1','1.5','2','3']]
 def half_beta(a,b):return mp.quad(lambda u:8*mp.power(mp.mpf('.5'),a)*mp.power(u,8*a-1)*mp.power(1-u**8/2,b-1),[0,1])
 for a in params:
  for b in params:
   integral=half_beta(a,b)+half_beta(b,a)
   close(integral,mp.gamma(a)*mp.gamma(b)/mp.gamma(a+b));betacases+=1
 result['beta_regularized_integral_cases']=betacases
 # Rational exponents: x=u^q removes the zero-end singularity entirely.
 reflectcases=0
 for q in range(2,10):
  for p in range(1,q):
   s=mp.mpf(p)/q
   integral=mp.quad(lambda u:q*u**(p-1)/(1+u**q),[0,1,mp.inf])
   close(integral,mp.pi/mp.sin(mp.pi*s));close(mp.gamma(s)*mp.gamma(1-s),integral);reflectcases+=1
 result['reflection_regularized_integral_cases']=reflectcases
 for z in [mp.mpc('.25','.3'),mp.mpc('.5','2'),mp.mpc('-1.3','.7'),mp.mpc('2','1'),mp.mpf('1.5')]:
  close(mp.gamma(z)*mp.gamma(1-z),mp.pi/mp.sin(mp.pi*z))
 gammacases=0
 for x in params+[mp.mpf(5),mp.mpf(8)]:
  integral=mp.quad(lambda u:4*u**(4*x-1)*mp.exp(-u**4),[0,1,2,mp.inf])
  close(integral,mp.gamma(x));close(mp.gamma(x+1),x*mp.gamma(x));gammacases+=1
  delta=mp.mpf('1e-6');R=max(mp.mpf(30),2*max(x-1,0));center=mp.quad(lambda t:t**(x-1)*mp.exp(-t),[delta,1,R])
  bound=delta**x/x+2*R**(x-1)*mp.exp(-R)
  require(0<=mp.gamma(x)-center<=bound)
 result['gamma_regularized_integral_cases']=gammacases
 bohrcases=0
 for j in range(1,11):
  x=mp.mpf(j)/10
  for n in [2,3,10,50,100,1000]:
   denom=mp.fprod(x+k for k in range(n));L=mp.factorial(n-1)*(n-1)**x/denom;U=mp.factorial(n-1)*mp.mpf(n)**x/denom;g=mp.gamma(x)
   require(L-mp.mpf('1e-65')<=g<=U+mp.mpf('1e-65'));close(U/L,(mp.mpf(n)/(n-1))**x);bohrcases+=1
 require(mp.e*mp.gamma(mp.mpf('1.25'))>1)
 result['bohr_mollerup_enclosures']=bohrcases
 inversions=[]
 for x in [mp.mpf('.1'),mp.mpf('.25'),mp.mpf(1),mp.mpf(4),mp.mpf(10)]:
  for T in [2,4,8]:
   approx=x**(-mp.mpf('.5'))*mp.quad(lambda t:mp.cos(t*mp.log(x))/mp.cosh(mp.pi*t),[0,T])
   err=abs(approx-1/(1+x));bound=2*x**(-mp.mpf('.5'))/mp.pi*mp.exp(-mp.pi*T)
   require(err<=bound);inversions.append({'x':str(x),'T':T,'error':mp.nstr(err,30),'analytic_tail_bound':mp.nstr(bound,30)})
 result['mellin_inverse_checks']=inversions
 I=mp.quad(lambda u:3/(1+u**3),[0,1,mp.inf]);close(I,2*mp.pi/mp.sqrt(3))
 U=mp.mpf(1000);cut=mp.quad(lambda u:3/(1+u**3),[0,1,10,100,U]);require(0<I-cut<3/(2*U*U))
 result['capstone_integral']={'exact':'2*pi/sqrt(3)','numeric':mp.nstr(I,50),'U':1000,'actual_cutoff_error':mp.nstr(I-cut,40),'analytic_tail_bound':'3/2000000'}
 # Principal value example uses an exact antiderivative at finite cutoffs.
 for c in [mp.mpf('.5'),mp.mpf(2),mp.mpf(3)]:
  R=mp.mpf('1e15');finite=mp.log((1+c*c*R*R)/(1+R*R))/2;require(abs(finite-mp.log(c))<mp.mpf('1e-28'))
 result['status']='PASS';result['assertions']=checks
 print(json.dumps(result,ensure_ascii=False,indent=2))
if __name__=='__main__':main()