#!/usr/bin/env python3
"""Recompute A01 asymptotic-integral examples. Python 3 + mpmath.
Run: python foundations-asymptotic-integrals-check.py [--output results.json]
High-precision quadrature is a diagnostic, not a certified interval routine.
The explicitly marked bounds follow from the algebraic inequalities in the pages.
"""
import argparse,json
from fractions import Fraction
import mpmath as mp
mp.mp.dps=70

def number(x):
    return mp.nstr(x,30)
def exact(x):
    return {'numerator':x.numerator,'denominator':x.denominator,'decimal':number(mp.mpf(x.numerator)/x.denominator)}
def check(condition,message):
    if not condition: raise AssertionError(message)

def run():
    report={'precision_decimal_digits':mp.mp.dps,'status':'passed','classification':{'numerical':'70-digit quadrature checks are evidence, not interval certificates','analytic':'Remainder bounds and rational endpoints are justified in the accompanying derivations'}}
    x=mp.mpf(8); fx=mp.quad(lambda t:mp.exp(-x*t)/(1+t),[0,1,mp.inf]); alternative=mp.exp(x)*mp.e1(x)
    check(abs(fx-alternative)<mp.mpf('1e-60'),'F integral versus exponential-integral identity')
    rows=[]
    for N in range(21):
        s=sum(((-1)**k*mp.factorial(k)/x**(k+1) for k in range(N)),mp.mpf(0))
        B=mp.factorial(N)/x**(N+1); r=fx-s
        check(abs(r)<=B,'F remainder bound N='+str(N));check((-1)**N*r>0,'F remainder sign N='+str(N))
        rows.append({'N':N,'sum':number(s),'signed_error':number(r),'absolute_error':number(abs(r)),'analytic_bound':number(B)})
    lo=sum((Fraction((-1)**k*__import__('math').factorial(k),8**(k+1)) for k in range(8)),Fraction())
    hi=sum((Fraction((-1)**k*__import__('math').factorial(k),8**(k+1)) for k in range(7)),Fraction())
    report['least_term']={'x':8,'value_quadrature':number(fx),'rows':rows,'exact_rational_enclosure':{'lower':exact(lo),'upper':exact(hi)}}
    # t=u^2 removes the integrable origin singularity before quadrature.
    x=mp.mpf(16);w=mp.quad(lambda u:2*mp.exp(-x*u*u)/(1+u*u),[0,1,mp.inf])
    w3=mp.sqrt(mp.pi)*(x**(-mp.mpf('.5'))-x**(-mp.mpf('1.5'))/2+3*x**(-mp.mpf('2.5'))/4)
    wb=mp.gamma(mp.mpf('3.5'))*x**(-mp.mpf('3.5'))
    check(0<w3-w<=wb,'Watson half-integer signed bound')
    report['watson']={'x':16,'value':number(w),'three_terms':number(w3),'absolute_error':number(abs(w-w3)),'analytic_bound':number(wb)}
    report['laplace']=[];report['stationary']=[]
    for L in [8,32,128]:
        l=mp.mpf(L);E=mp.quad(lambda t:mp.exp(-l*(t+t*t)),[0,1]);E2=1/l-2/l**2
        J=mp.quad(lambda t:mp.exp(-l*(t*t/2+t**4/4)),[-1,0,1]);J2=mp.sqrt(2*mp.pi/l)*(1-3/(4*l))
        JB=mp.sqrt(2*mp.pi/l)*105/(32*l*l)+2*mp.exp(-l/2)/l
        check(abs(J-J2)<=JB,'interior Laplace two-term analytic bound')
        report['laplace'].append({'lambda':L,'endpoint_exact_numeric':number(E),'endpoint_two_terms':number(E2),'endpoint_error':number(abs(E-E2)),'interior_exact_numeric':number(J),'interior_two_terms':number(J2),'interior_error':number(abs(J-J2)),'interior_analytic_bound':number(JB)})
        K=mp.sqrt(mp.pi)/mp.sqrt(1-1j*l/2);K1=mp.sqrt(2*mp.pi/l)*mp.exp(1j*mp.pi/4)
        # Oscillatory quadrature with phase-aware short panels on [-12,12].
        # Omitted Gaussian tail <= exp(-144)/12, below 1e-60.
        panels=[mp.mpf(j)/8 for j in range(-96,97)]
        Kq=mp.quad(lambda t:mp.exp(-t*t+1j*l*t*t/2),panels)
        check(abs(Kq-K)<mp.mpf('1e-58'),'stationary exact complex formula versus quadrature')
        report['stationary'].append({'lambda':L,'exact_real':number(K.real),'exact_imag':number(K.imag),'quadrature_difference':number(abs(Kq-K)),'leading_relative_error':number(abs(K-K1)/abs(K))})
    # Nonstationary endpoint formula with a nontrivial amplitude, independent quadrature.
    l=mp.mpf(7);Q=mp.quad(lambda t:t*mp.exp(1j*l*t),[0,1]);Qexact=mp.exp(1j*l)/(1j*l)+(mp.exp(1j*l)-1)/l**2
    check(abs(Q-Qexact)<mp.mpf('1e-60'),'nonstationary two-term exact formula')
    report['nonstationary']={'lambda':7,'exact_formula_difference':number(abs(Q-Qexact))}
    # Exact quartic scaling on R: numerical integral after the analytic rescaling.
    C=mp.quad(lambda u:mp.exp(-u**4),[-mp.inf,0,mp.inf]);Cexact=mp.gamma(mp.mpf('.25'))/2
    check(abs(C-Cexact)<mp.mpf('1e-60'),'quartic concentration constant')
    report['degenerate_quartic']={'constant_numeric':number(C),'constant_gamma':number(Cexact),'scale':'lambda^(-1/4), not lambda^(-1/2)'}
    return report
if __name__=='__main__':
    parser=argparse.ArgumentParser();parser.add_argument('--output');args=parser.parse_args()
    payload=json.dumps(run(),ensure_ascii=False,indent=2)
    if args.output:
        with open(args.output,'w') as handle:handle.write(payload+'\n')
    else:print(payload)