Skip to main content

bessel_in_pair_impl

Function bessel_in_pair_impl 

Source
pub fn bessel_in_pair_impl<P, E, V, const NS: usize, const NL: usize, const NF: usize, const N: i32, const SCALED: bool>(
    x: V,
    t: &BesselI<E, NS, NL, NF>,
) -> (V, V)
Expand description

I_N(x) for N >= 2, or e^{-|x|} I_N(x) when SCALED, by downward recurrence on the ratios rather than on the values.

§Why ratios

Writing r_k = I_k(x)/I_{k-1}(x), the three-term recurrence I_{k-1} = I_{k+1} + (2k/x) I_k divides through to

r_k = \frac{1}{2k/x + r_{k+1}}

which is the same continued fraction Boost evaluates with Lentz’s method in CF1_ik. The point for a vector unit is that every r_k lies in (0, 1), so nothing can overflow and no rescaling is needed anywhere. The textbook alternative, Miller’s linear downward pass on the values themselves, has intermediates growing like 2^M M!/x^M, which leaves f64 range around order 50 and f32 range around order 8, and needs a per-lane rescale select inside the loop to survive. That is three selects per iteration to buy nothing.

Seeding is from I_0, which the closed form already provides, and I_N = I_0 \prod r_k. The scaled and unscaled forms differ only in that seed. A ratio is scale-free, so SCALED never reaches the loop.

§Why forward recurrence is not used

The textbook rule is “forward when N < x, downward otherwise”, which would bound both trip counts by N alone. Measured, it does not work: forward recurrence’s amplification grows with N faster than it decays with x, giving 2.1e-06 at N = 50, x = 100, and N = 80 never reaches 1e-13 for any x up to 700. No crossover rescues it: the best fitted rule still left 1.4e-03. So there is one path here, and therefore no select between paths at all.

§Trip count

N + 24 + coeff * x, with the per-lane start following that lane’s own x. Lanes are free to start higher than they need, because the recurrence is self-correcting downward from a zero seed. No lane is ever cut short, and the loop simply runs until every lane has walked down to k = 1. The packet therefore pays its worst lane, which is the standing trade for a data-dependent trip count here.

Last built: 2026-09-08 21:35:55 UTC