Finite exact-interval evaluation of scaled spherical Bessel values

Inputs are those of Scaled spherical Bessel evaluation on positive real arguments, a nonempty finite list of valid targets, finite exact positive rational absolute tolerances, optional exact positive rational relative tolerances, and finite integer caps $W,Q,M_{\max },p_{\max }\ge 1$. $W$ is the maximum exact-arithmetic work count, $Q$ the maximum stage and represented-point interface index, $M_{\max }$ the maximum Taylor degree, and $p_{\max }$ the maximum final dyadic precision. A Gaussian rational is $a+ib$ with $a,b\in \mathbb{Q}$. A Laurent polynomial is a finite sum ${\sum }_k{a}_k{z}^k$ with integer exponents and Gaussian-rational coefficients; its formal derivative uses $(z^k{)}^{\mathrm{\prime }}=k{z}^{k-1}$. An expression graph below is a finite directed acyclic graph of exact rational arithmetic operations.

Factorials are those of The factorial $n!$ and the falling factorial $n^{\underset{‾}{k}}$, defined by recursion in $\mathbb{N}$; they are nonzero by that definition and Canonical naturals are positive and strictly increasing, so every displayed factorial denominator is defined. For a rational $x$, $\lceil x\rceil$ denotes the least integer greater than or equal to $x$, computed by exact integer division.

Reject an invalid order, target, tolerance, cap, or interval input not contained in $(0,\mathrm{\infty })$ with domain. For a represented point, an interval returned by the input interface may meet zero; such a stage is skipped and a later interface index is tried.

Construct Laurent polynomials by $A_0=0,B_0=z^{-1},C_0=-i{z}^{-1},$ $A_{r+1}=r{z}^{-1}{A}_r-A_r^{\mathrm{\prime }}-B_r,B_{r+1}=r{z}^{-1}{B}_r-B_r^{\mathrm{\prime }}+A_r,C_{r+1}=r{z}^{-1}{C}_r-C_r^{\mathrm{\prime }}-i{C}_r$ for $r=0,\dots ,\mathrm{\ell }-1$. With $s$ and $c$ denoting the sine and cosine input enclosures, form the unscaled outputs $u=A_{\mathrm{\ell }}c+B_{\mathrm{\ell }}s,v=(A_{\mathrm{\ell }}^{\mathrm{\prime }}+B_{\mathrm{\ell }})c+(B_{\mathrm{\ell }}^{\mathrm{\prime }}-A_{\mathrm{\ell }})s,w=C_{\mathrm{\ell }}(c+is),t=(C_{\mathrm{\ell }}^{\mathrm{\prime }}+i{C}_{\mathrm{\ell }})(c+is)\mathrm{.}$ Form the scaled outputs $c_{\mathrm{\ell }}{z}^{-\mathrm{\ell }}u,c_{\mathrm{\ell }}{z}^{-\mathrm{\ell }}v,d_{\mathrm{\ell }}^{-1}{z}^{\mathrm{\ell }+1}{C}_{\mathrm{\ell }},d_{\mathrm{\ell }}^{-1}{z}^{\mathrm{\ell }+1}(C_{\mathrm{\ell }}^{\mathrm{\prime }}+i{C}_{\mathrm{\ell }})\mathrm{.}$ Collect like Laurent powers exactly; in particular cancel negative powers in the scaled $h_{\mathrm{\ell }}^{(1)}$ value when they cancel algebraically. A remaining $z^{-1}$ in the scaled derivative is evaluated only on positive input intervals.

For real intervals use exact rational endpoint hulls: $[a,b]+[c,d]=[a+c,b+d]$, $[a,b][c,d]=[\min (ac,ad,bc,bd),\max (ac,ad,bc,bd)]$, and $[a,b{]}^{-1}=[1\mathrm{/}b,1\mathrm{/}a]$ when $0<a\le b$. For complex rectangles, apply these real interval operations to the real and imaginary components of $(x+iy)(u+iv)=(xu-yv)+i(xv+yu)$. Thus interval arithmetic encloses every pointwise sum and product; division is only by a positive power of the positive argument or by a positive integer constant. Compute powers, factorials, and rational endpoints by finite exact arithmetic.

At a positive current box $I=[a,b]$, choose $M\ge 0$ with $M+2\ge 2b$ and form the rational interval polynomial and rational tail $T_M(I)={\sum }_{m=0}^M\frac{(iI{)}^m}{m!},E_M=\frac{2b^{M+1}}{(M+1)!}\mathrm{.}$ The interval power is repeated complex-rectangle multiplication. Enclose $\cos (I)$ by $\mathrm{Re}{T}_M(I)+[-E_M,E_M]$ and $\sin (I)$ by $\mathrm{Im}{T}_M(I)+[-E_M,E_M]$. Evaluate each requested output graph with these leaves and $z\in I$, and round each final endpoint outward to the grid $2^{-p}\mathbb{Z}$.

Use the following deterministic finite search. Charge the Laurent recurrence, graph construction, and coefficient collection to the work counter as they are performed. For stages $q=1,\dots ,Q$, first take the fixed interval unchanged, or call the represented-point interface at index $q$. Skip a represented-point stage unless its returned interval has $a>0$. Set $p=q$ and $M=\max (q,\lceil 2b\rceil )$. If $p>p_{\max }$ or $M>M_{\max }$, stop with inconclusive. Evaluate all requested graphs and the displayed $T_M,E_M$ using the work counter: each exact rational or Gaussian-rational arithmetic operation or comparison, each Taylor term, each endpoint-rounding operation, and each interface call costs one unit. If the next operation would exceed $W$, stop with inconclusive before performing it. If the output format cannot represent a certified endpoint, return range with no finite claim. If every requested target has radius at most its absolute tolerance, return success; for each optional relative tolerance, set its flag true exactly when the enclosure excludes zero and the relative test in the contract passes. A failed relative test does not prevent absolute success. If no stage succeeds before a cap or the search ends, return inconclusive. No floating-point recurrence is used.