Module ratio 

A_nu(x) = I_nu(x) / I_{nu-1}(x), the modified Bessel ratio, and its inverse.

Where it comes from

With p = 2 nu this is the mean resultant length of a von Mises-Fisher distribution on the sphere S^{p-1} as a function of its concentration kappa, and its inverse is the maximum-likelihood concentration from an observed mean resultant length: the one step of every vMF fit that is not a matrix product. p = 2 (nu = 1) is the von Mises circle, I_1/I_0. p = 3 (nu = 3/2) collapses to the elementary Langevin function coth x - 1/x, which is the langevin family. The order arrives as a plain vector, so those cases reach the integer and half-integer Bessel kernels through the order simplifier, not the general real-order machinery.

Forward, in three arms

The obvious quotient of the two scaled Bessel functions fails at small x when the order is large: e^{-x} I_nu(x) is (x/2)^nu / Gamma(nu+1) there and underflows to zero (I_150(0.5) is 1e-457) while the ratio, about x / 2nu, is ordinary. And it keeps failing well past that: the scaled I_nu(x) is below e^{-0.17 nu} at x = nu.

• x <= 0.9 sqrt(nu): the ratio of the two power series, A = (x / 2nu) S_nu / S_{nu-1} with S_a = sum (x^2/4)^k / (k! (a+1)_k), twelve terms each (1 / (4^12 12!) at the edge). Never underflows: only the ratio is formed.
• x < 8 nu: the Perron continued fraction for I_{nu}/I_{nu-1} (cf1_i_ratio, the Bessel kernels’ own), which converges for every (nu, x) in about 6 sqrt(x) iterations.
• else: the quotient of the scaled Bessels from [bessel_iv], where x dominates the order and both are ordinary numbers, and the asymptotic arm inside is fast.

The derivative

From the recurrence I_{nu-1} - I_{nu+1} = (2 nu / x) I_nu: A' = 1 - A^2 - (2 nu - 1) A / x, a closed form in A itself, which is what makes the inverse a one-evaluation Newton and gives Dual its factor for free.

Inverse

Banerjee’s kappa_0 = r (p - r^2) / (1 - r^2) (within ten percent everywhere), then newtons_method on A(kappa) - r with the derivative above, bracketed by a factor of two, eight iterations at most. Sra (2012) found two steps from that seed reach working precision, and the tolerance is a few ulp of r. Below r = 1e-8 the answer is 2 nu r outright (the next term is r^3).

The inverse is ill-conditioned as r -> 1, where kappa ~ (p-1)/(2(1-r)): an ulp of r moves kappa by 2 kappa^2 eps / (p-1), so the relative error grows like kappa, and the result is the exact inverse of the given r only to that extent. A complement form taking 1 - r, the way inv_langevin_1m does for p = 3, is the remedy and is not built.

Functions

bessel_i_ratio_1m_impl

1 - A_nu(x), the complement of the ratio, accurate where A is within an ulp of 1.

bessel_i_ratio_deriv

A' = 1 - A^2 - (2 nu - 1) A / x from A itself.

bessel_i_ratio_deriv_1m

A' from the complement c = 1 - A, with 1 - A^2 as c (2 - c): no cancellation where A is within an ulp of 1.

bessel_i_ratio_impl

I_nu(x) / I_{nu-1}(x), for x >= 0 and nu >= 1. Odd in x.

inv_bessel_i_ratio_1m_impl

The kappa with 1 - A_nu(kappa) = t, for 0 < t <= 2 (t = 1 - r), the complement form of inv_bessel_i_ratio_impl: well conditioned as t -> 0, where the plain form loses 2 kappa eps / (p - 1) to the rounding of r.

inv_bessel_i_ratio_impl

The kappa >= 0 with I_nu(kappa) / I_{nu-1}(kappa) = r, for 0 <= r < 1. Odd in r.