Modified Bessel function of the first kind, order 0: $I_0(x)$.

Composite implementation via the Abramowitz & Stegun polynomial approximation (Table 9.8.1):

• |x| <= 3.75 — polynomial series in $(x / 3.75)^2$.
• |x| > 3.75 — polynomial series in $3.75 / |x|$, multiplied by $e^{|x|} / \sqrt{|x|}$ to recover the exponential growth.

Each branch runs on input clamped into its own interval, so the one where does not pick contributes no nan to the gradient.

A float64 input takes the Cephes double-precision series instead.