The functions \( I_0(x) \) and \( I_1(x) \) solve the modified Bessel equation
\[ x^2 \frac{d^2 y}{dx^2} + x \frac{dy}{dx} - (x^2 + \nu^2)\, y = 0 \]
and grow exponentially with \( x \). They are obtained from the ordinary Bessel functions through \( I_\nu(x) = i^{-\nu} J_\nu(i x) \), and admit the series
\[ I_n(x) = \sum_{m=0}^{\infty} \frac{1}{m!\,(m+n)!} \left(\frac{x}{2}\right)^{2m+n}, \]
with the large- \( x \) asymptotic \( I_\nu(x) \sim e^{x} / \sqrt{2\pi x} \).
Function shape
Domain and special values
| \( x \) | \( I_0(x) \) | \( I_1(x) \) |
|---|---|---|
| \( 0 \) | \( 1 \) | \( 0 \) |
| \( -x \) | \( I_0(x) \) (even) | \( -I_1(x) \) (odd) |
| \( +\infty \) | \( +\infty \) | \( +\infty \) |
Both functions are defined for all real \( x \) (the implementation works on \( |x| \), restoring the sign for the odd \( I_1 \)).
Usage
Accuracy and performance
besseli0/besseli1 switch at \( x = 7.75 \): a minimax polynomial in \( (x/2)^2 \) below, and a scaled \( \sqrt{x}\,e^{-x} I_\nu(x) \) minimax form (Remez-fitted) above. Compared with the netlib RIBESL reference they are dramatically faster (hundreds of ns/eval → ~10 ns/eval) at matching accuracy:
