Finite certified evaluation of integer-order Bessel functions

Inputs are an integer order $n\ge 0$, a target $J$ or scaled $I$, an exact dyadic real argument $x$, positive rational absolute tolerance ${\epsilon }_a$, optional positive rational relative tolerance ${\epsilon }_r$, a fixed exponent range $(e_{\min },e_{\max })$, and finite positive integer caps $S_{\max },K_{\max },M_{\max },N_{\max },p_{\max },W$, with $p_{\max }\ge 2$. The caps bound stages, series terms, midpoint panels, Miller terminal index, dyadic precision, and exact interval operations and subroutine calls. Run stages $q=1,\dots ,\min (S_{\max },K_{\max },M_{\max },p_{\max })$, using precision $p=\max (2,\min (q,p_{\max }))$ and elementary-range truncation $q$. Reject $n<0$, invalid targets, nonpositive tolerances, invalid caps, $e_{\min }>e_{\max }$, or $x<0$ for the scaled $I$ target with domain. The scaled target is always $e^{-x}{I}_n(x)$.

Use exact rational endpoint formulas for finite sums, recurrences, mesh construction, and error bounds, with final endpoints rounded outward to precision $p$. Charge every exact arithmetic operation, range-evaluator node operation, term, rounding, and loop update to the shared counter $W$; interrupt a nested elementary-range call with inconclusive before its next operation would exceed the remaining budget. Factorials use The factorial $n!$ and the falling factorial $n^{\underset{‾}{k}}$, defined by recursion in $\mathbb{N}$. Call Certified range evaluation for elementary real and complex functions for the real sine, cosine, and exponential ranges in the midpoint branches. Before each call, round its rational input interval outward to the fixed dyadic format. Its finite outputs are sound by Certified real elementary-function evaluations enclose their ranges. Count each exact operation, term, rounding, and subroutine call against $W$; if the next operation would exceed $W$, return inconclusive. A non-finite elementary call makes that branch unavailable. Final endpoint overflow returns range.

Before forming any stage constant, if the target is scaled $I$ and $x=0$, use precision $p=p_{\max }$, form $[1,1]$ when $n=0$ and $[0,0]$ when $n\ge 1$, and round its endpoints outward to that format. If either conversion overflows, return range. Otherwise apply the ordinary absolute-radius test to the rounded candidate and return success only when it passes; return inconclusive if it fails. The internal exact candidate has radius zero, but outward underflow may widen its stored enclosure. Set a successful output's relative flag by the usual contract test, so a zero value receives no relative certificate. This branch does not call the elementary-range procedure.

At stage $q$, compute an outward dyadic interval $P_q$ containing $\pi$ from the Machin identity in Trigonometric identities, Machin's formula and reduced Taylor tails and the arctangent series and geometric-tail argument in Range-reduced exponential, logarithm and arctangent series have certified tails: for $r=1\mathrm{/}5,1\mathrm{/}239$, use the $q$-term sum ${\sum }_{k=0}^{q-1}(-1{)}^k{r}^{2k+1}\mathrm{/}(2k+1)$ and tail bound $r^{2q+1}\mathrm{/}((2q+1)(1-r^2))$, combine by $16\arctan (1\mathrm{/}5)-4\arctan (1\mathrm{/}239)$, and round outward to the selected format at precision $p$. If this conversion is non-finite or $\inf P_q\le 0$, midpoint candidates are unavailable at this stage; the independent power-series or Miller branches may still be attempted.

For a requested $J_n(x)$, first form a power-series candidate if possible. For any order $r\ge 0$, define the finite subroutine $\mathrm{SeriesJ}(r)$: start $a_0=(x\mathrm{/}2{)}^r\mathrm{/}(0!r!)$ and $a_{k+1}=-a_k{x}^2\mathrm{/}(4(k+1)(k+r+1))$; sum $a_0,\dots ,a_{q-1}$, and compute $a_q$. Put ${\rho }_J=\frac{x^2}{4(q+1)(q+r+1)}\mathrm{.}$ If ${\rho }_J<1$, return the interval partial sum enlarged by $\mathrm{\mid }{a}_q\mathrm{\mid }\mathrm{/}(1-{\rho }_J)$; otherwise mark this series candidate unavailable. The tail lemma Convergence and geometric tails for integer-order Bessel series proves each returned interval contains $J_r(x)$.

For each requested target, also form its composite-midpoint candidate with $m=q$ panels, $h=\pi \mathrm{/}m$, and nodes ${\theta }_j=(j+\frac{1}{2})h$, $j=0,\dots ,m-1$, using $P_q$ as the interval for $\pi$. For $J_n$, call the elementary-range algorithm for $\sin ({\theta }_j)$, form $n{\theta }_j-x\sin ({\theta }_j)$, and call it for the cosine of that interval. Let $Q_m$ be the interval sum $h{\sum }_j\cos (n{\theta }_j-x\sin {\theta }_j)$ and $E_J=\frac{h_{\max }^{}}{24}((n+\mathrm{\mid }x\mathrm{\mid }{)}^2+\mathrm{\mid }x\mathrm{\mid }),h_{\max }=\sup (P_q)\mathrm{/}m\mathrm{.}$ The candidate is $Q_m\mathrm{/}{P}_q+[-E_J,E_J]$. For the scaled $I$ target with $x>0$, call the elementary algorithm for $\cos ({\theta }_j)$, form $-x(1-\cos ({\theta }_j))$, call it for the real exponential, and also call it for $\cos (n{\theta }_j)$. Let $Q_m$ be $h$ times the sum of the resulting product intervals and set $E_I=\frac{h_{\max }^{}}{24}((x+n{)}^2+x)\mathrm{.}$ The candidate is $Q_m\mathrm{/}{P}_q+[-E_I,E_I]$. The exact midpoint error lemmas Composite midpoint enclosure for the integer Bessel integral and Composite midpoint enclosure for the scaled modified Bessel integral prove these inclusions.

For a $J_n(x)$ Miller candidate, only try this branch when $x>0$ and $N=n+q\le N_{\max }$. Obtain $\mathrm{SeriesJ}(r)$ for $r=0,1,n,N+1$; if any is unavailable, skip the branch. Set $D_0=-J_1$, $D_1=J_0$, and compute $D_{j+1}=2j{D}_j\mathrm{/}x-D_{j-1}$ for $j=1,\dots ,N$ by interval arithmetic. Set $u_{N+1}=0,u_N=1$, and compute $u_{j-1}=2j{u}_j\mathrm{/}x-u_{j+1}$ for $j=N,\dots ,1$. Require the intervals for $D_{N+1}$, $u_s$, $J_s$, and $J_s+q_M{D}_s$ below to exclude zero. Form $q_M=-J_{N+1}\mathrm{/}{D}_{N+1}$. In the fixed order $s=0,1$, whenever the tests pass, form ${\widehat{J}}_n=\frac{J_s{u}_n}{u_s},E_M=\frac{q_M(J_s{D}_n-J_n{D}_s)}{J_s+q_M{D}_s},C_M={\widehat{J}}_n-E_M\mathrm{.}$ The exact finite terminal-contamination identity Exact finite terminal-contamination error identity proves that each $C_M$ contains $J_n(x)$. The Casoratian lemma Positive discrete Wronskian for the Bessel recurrence pair supplies the independent recurrence pair and positive invariant required by that identity.

Collect every available sound candidate for the requested target. If more than one is available, intersect them; if the intersection is empty, discard this stage. Round the selected candidate outward to precision $p$. Return success at the first stage whose radius is at most ${\epsilon }_a$. Set relative_certified to true exactly when relative tolerance was requested, zero is excluded, and the contract's relative-radius test passes; failure of that relative test does not prevent absolute success. If no stage succeeds before the finite caps, return inconclusive. No branch uses a limit in the Miller terminal index.