Soundness and eventual absolute precision for spherical Bessel enclosures

For $y_{\mathrm{\ell }}=(-z{)}^{\mathrm{\ell }}{D}^{\mathrm{\ell }}f$, differentiation as in [F1] gives $y_{\mathrm{\ell }+1}=\mathrm{\ell }{y}_{\mathrm{\ell }}\mathrm{/}z-y_{\mathrm{\ell }}^{\mathrm{\prime }}$. Starting from $j_0=\sin z\mathrm{/}z=A_0\cos z+B_0\sin z$ with $A_0=0,B_0=z^{-1}$, the derivative formula $(A_{\mathrm{\ell }}\cos z+B_{\mathrm{\ell }}\sin z{)}^{\mathrm{\prime }}=(A_{\mathrm{\ell }}^{\mathrm{\prime }}+B_{\mathrm{\ell }})\cos z+(B_{\mathrm{\ell }}^{\mathrm{\prime }}-A_{\mathrm{\ell }})\sin z$ gives exactly the algorithm's $A_{\mathrm{\ell }+1},B_{\mathrm{\ell }+1}$. Starting from $h_0^{(1)}=-i{e}^{iz}\mathrm{/}z=C_0{e}^{iz}$ with $C_0=-i\mathrm{/}z$, the componentwise product rule [F5]–[F6] gives $(C_{\mathrm{\ell }}{e}^{iz}{)}^{\mathrm{\prime }}=(C_{\mathrm{\ell }}^{\mathrm{\prime }}+i{C}_{\mathrm{\ell }})e^{iz}$ and exactly the algorithm's recurrence for $C_{\mathrm{\ell }+1}$. Euler's identity [F6] identifies the algorithm's $c+is$ with $e^{iz}$. Thus $u,v,w,t$ represent $j_{\mathrm{\ell }},j_{\mathrm{\ell }}^{\mathrm{\prime }},h_{\mathrm{\ell }}^{(1)},(h_{\mathrm{\ell }}^{(1)}{)}^{\mathrm{\prime }}$, respectively.

The two series for $e^{-iz}$ and $e^{iz}$ are absolutely convergent by [F2], so [F7] identifies their product with the sum of the Cauchy-product coefficients. For total degree $n$, its coefficient is $z^n{S}_n$, where $S_n={\sum }_{r=0}^n(-i{)}^r{i}^{n-r}\mathrm{/}(r!(n-r)!)$. For $1\le r\le n$, combining over the denominator $r!(n+1-r)!$ and using $(n+1-r)!=(n+1-r)(n-r)!$ and $r!=r(r-1)!$ gives $1\mathrm{/}(r!(n+1-r)!)=(1\mathrm{/}(n+1))(1\mathrm{/}(r!(n-r)!)+1\mathrm{/}((r-1)!(n+1-r)!))$. For $r=0$ or $r=n+1$, omit the summand whose index is outside $0,\dots ,n$; the remaining summand equals the left side by $0!=1$. Splitting the finite sum for $S_{n+1}$ with this identity and reindexing its second part gives $S_{n+1}=(i{S}_n-i{S}_n)\mathrm{/}(n+1)=0$ for every $n\ge 0$, while $S_0=1$. The factorials and $n+1$ are nonzero by [F2], so these divisions are defined. Hence the product series has constant coefficient $1$ and every positive-degree coefficient $0$, so $e^{-iz}{e}^{iz}=1$.

Fix a positive input box $I=[a,b]$ and choose $M$ with $M+2\ge 2b$. For every $z\in I$, the first omitted absolute exponential-series term is at most $b^{M+1}\mathrm{/}(M+1)!$. Each subsequent ratio is at most $b\mathrm{/}(M+2)\le 1\mathrm{/}2$, so the sum of the omitted terms is at most $2b^{M+1}\mathrm{/}(M+1)!=E_M$. Thus the real and imaginary parts of $T_M(z)$, enlarged by $[-E_M,E_M]$, contain $\cos z$ and $\sin z$. Exact rational interval addition and multiplication enclose pointwise addition and multiplication by their endpoint-hull formulas; inversion of a positive interval $[a,b]$ is $[1\mathrm{/}b,1\mathrm{/}a]$. Evaluation over the finite graph and outward rounding therefore preserve containment of every graph value on $I$.

Each cap bounds its corresponding loop or finite expansion. The Laurent graph for fixed $\mathrm{\ell }$ is finite; each rational operation and each terminating represented-point interface call takes finite work, and every cap exhaustion has an explicit status return. Therefore every capped execution terminates, including one that never obtains a positive point box or never passes its requested tests.

Fix a positive represented point $z_0$. By [F10], choose a sufficiently large index $q_0$ that $2^{-q_0}<z_0\mathrm{/}4$. Write $I_{q_0}=[a_0,b_0]$; the interface width bound gives $a_0>0$. For all $q\ge q_0$, $I_q\subset [a_0\mathrm{/}2,b_0+1]=:[a,B]\subset (0,\mathrm{\infty })$, because $I_q$ contains $z_0$, has width at most $2^{-q_0}$, and $z_0-2^{-q_0}>a_0\mathrm{/}2$ while $z_0+2^{-q_0}<b_0+1$. Thus $M=q$ eventually and $E_M\le 2B^{q+1}\mathrm{/}(q+1)!\to 0$, since the successive ratios are eventually at most $1\mathrm{/}2$. Put ${\delta }_q=\mathrm{width}(I_q)$. For $m\ge 1$, positivity of $I_q=[a_q,b_q]$ makes repeated complex-rectangle multiplication of $i{I}_q$ exact along one coordinate: its $m$th power is $i^m[a_q^m,b_q^m]$ up to the sign absorbed by the interval endpoints. Therefore $\mathrm{width}(I_q^m)=b_q^m-a_q^m\le m{B}^{m-1}{\delta }_q$ by factoring the difference of powers. Each degree-$m$ term in $T_M(I_q)$ has component width at most ${\delta }_q{B}^{m-1}\mathrm{/}(m-1)!$, so each component of the polynomial interval has width at most ${\delta }_q{\sum }_{m=1}^M{B}^{m-1}\mathrm{/}(m-1)!$. The terms of the last series are nonzero because $B>0$, and their successive ratio is $B\mathrm{/}m\to 0$; [F8] makes its sum finite, say $C_B$. With the Taylor tail added, the sine and cosine enclosure widths tend to zero.

Multiplying the unscaled expressions of step 1.1 by the scales in Scaled spherical Bessel evaluation on positive real arguments gives the four scaled identities. For $H_{\mathrm{\ell }}^{\mathrm{s}\mathrm{c}}$ and $D{H}_{\mathrm{\ell }}^{\mathrm{s}\mathrm{c}}$, cancel $e^{-iz}{e}^{iz}$ by step 1.2; the derivative scale acts on the original $h_{\mathrm{\ell }}^{(1)}$ derivative as defined, not on $H_{\mathrm{\ell }}^{\mathrm{s}\mathrm{c}}$. Starting with $C_0=-i{z}^{-1}$, if $C_{\mathrm{\ell }}$ has exponents $-1,\dots ,-\mathrm{\ell }-1$, then $\mathrm{\ell }{z}^{-1}{C}_{\mathrm{\ell }}$ and $C_{\mathrm{\ell }}^{\mathrm{\prime }}$ have exponents at least $-\mathrm{\ell }-2$, while $-i{C}_{\mathrm{\ell }}$ has exponents at least $-\mathrm{\ell }-1$; induction gives only exponents $-1,\dots ,-\mathrm{\ell }-1$ in $C_{\mathrm{\ell }}$. Thus $z^{\mathrm{\ell }+1}{C}_{\mathrm{\ell }}$ is a polynomial and $z^{\mathrm{\ell }+1}(C_{\mathrm{\ell }}^{\mathrm{\prime }}+i{C}_{\mathrm{\ell }})$ has at worst exponent $-1$, which is well-defined on the positive input domain.

The remaining expression graph is fixed and finite. On $[a,B]$, each reciprocal has width ${\delta }_q\mathrm{/}(a_q{b}_q)\le {\delta }_q\mathrm{/}{a}^2$ and is bounded. Sums of bounded intervals with widths tending to zero also have widths tending to zero; for bounded real intervals $X,Y$, $\mathrm{width}(XY)\le \sup \mathrm{\mid }X\mathrm{\mid }\mathrm{width}(Y)+\sup \mathrm{\mid }Y\mathrm{\mid }\mathrm{width}(X)$. Applying these bounds to the real component operations in each complex rectangle proves by induction that every output rectangle or interval has width tending to zero. By [F10], the outward rounding increment $2^{1-q}$ also tends to zero, so all output radii tend to zero. Thus every positive rational absolute tolerance is passed at some finite stage, and the finite work through that stage gives finite sufficient caps. Constant interfaces prove the exact positive-dyadic case as well.

By steps 1.1 and 2.1 the graph roots are the exact requested targets. Step 1.3 encloses the sine and cosine leaves for every admissible argument; induction through the finite graph, using the interval formulas in [F4], encloses every sum, product, reciprocal, and requested output value. Outward endpoint rounding preserves these containments. The success-radius and relative-flag tests are precisely the contract in [F9], so every finite success is sound for point or whole-interval input.

If a requested target has value $v\mathrm{\ne }0$, any sound enclosure $B$ of step 3.1 has diameter $2\mathrm{rad}(B)$ and contains $v$, so $d(B,0)\ge \mathrm{\mid }v\mathrm{\mid }-2\mathrm{rad}(B)$. Eventually $d(B,0)\ge \mathrm{\mid }v\mathrm{\mid }\mathrm{/}2>0$ and $\mathrm{rad}(B)\mathrm{/}d(B,0)\to 0$, giving any requested positive rational relative tolerance along the refined sequence. The algorithm may already have stopped with absolute success before that relative test passes. To obtain the relative flag, rerun with successively smaller positive rational absolute tolerances and increasing finite caps: step 2.2 supplies a successful finite capped run for each such tolerance, and its returned radius is bounded by that tolerance. For a finite list of nonzero requested targets, sufficiently small absolute tolerances make all requested relative flags true. If $v=0$, every sound enclosure contains zero and no relative flag is available. A fixed interval input encloses its whole image, which may have positive diameter, so the point-mode limiting argument asserts no arbitrary-radius guarantee there.