Module ik 

Modified Bessel functions of the first kind, orders 0 and 1, scaled and unscaled.

Two regions, not four

$I_0$ and $I_1$ split at x = 7.75 and nowhere else: an ascending series in $a = x^2/4$ below, and $e^{x} P(1/x)/\sqrt{x}$ above. That is the whole shape.

The obvious reference for a SIMD Bessel is fdlibm, and fdlibm is the wrong model. Its j0f splits the asymptotic envelope alone into four sub-intervals with a rational apiece, which is optimal when a branch picks one and skips the rest, and pathological here. A vector unit evaluates all four and discards three. The tables underneath this kernel take Boost’s route instead: fewer regions, higher degree in each. Same trade the Faddeeva kernel makes for the same reason.

The scalar sources carry a third region near the top of the range. It is not an accuracy region. It exists so $e^x$ cannot overflow before $/\sqrt{x}$ brings the product back down, and is reached here only under [thermite::unlikely], so the common path pays one compare.

f32 $I_0$ is the one exception, and there the far fit is genuine: its large minimax is fitted over [7.75, 50] and its constant term is wrong in the seventh digit, so extending it to infinity costs about 16 ulp. See [crate::tables::bessel].

The scaled forms are the cheaper ones

$e^{-x} I_n(x)$ is not a wrapper that multiplies an exponential back out. Above 7.75 the tables are the scaled value, so the scaled entry points skip the exponential entirely and the unscaled ones pay for it. Below 7.75 the relationship inverts. Each form therefore costs one transcendental in exactly one of its two arms, and neither is built from the other. That is the relationship zetac has to zeta, for the same reason ($I_0(800)$ overflows f64 while $e^{-800}I_0(800)$ is a perfectly ordinary 0.0141).

Functions

asymptotic_series_g

asymptotic_series_v over a general arithmetic: the argument z in C, the order in its real primal R, and the far mask supplied by the caller.

bessel_i0_impl

$I_0(x)$, or $e^{-|x|} I_0(x)$ when SCALED.

bessel_i1_impl

$I_1(x)$, or $e^{-|x|} I_1(x)$ when SCALED.

bessel_in_impl

$I_N(x)$ for N >= 2, or $e^{-|x|} I_N(x)$ when SCALED.

bessel_in_pair_impl

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.

bessel_iv_impl

I_n(x) with a per-lane order, or e^{-|x|} I_n(x) when SCALED.

bessel_k0_impl

$K_0(x)$, or $e^{x} K_0(x)$ when SCALED.

bessel_k1_impl

$K_1(x)$, or $e^{x} K_1(x)$ when SCALED.

bessel_kn_recur

$K_N(x)$ for N >= 2 from $K_0$ and $K_1$, by upward recurrence.

bessel_kv_impl

K_n(x) with a per-lane order, or e^{x} K_n(x) when SCALED.

unscale_i_pair_masked

unscale_i_pair with the far mask supplied, for an arithmetic where “far” is not a plain comparison: over C it is Re z >= threshold, and z itself is the exponent.