Module langevin 

Langevin function L(x) = coth(x) - 1/x and its inverse, shared by every real element type. The per-precision pieces (the polynomial tables, and Newton vs Halley for the inverse) come in from ps.rs/pd.rs.

Forward

coth(x) and 1/x both grow like 1/x while their difference is only x/3, so the direct form loses relative accuracy as 3u/x^2. That is not a corner case, at x = 0.1 it is already 300 ulp. Below the crossover X0 = 2 the function is therefore an odd polynomial x * p(x^2) (minimax, fitted to L(x)/x), and above it

L(x) = 1 - \frac{1}{x} + \frac{2q}{1 - q}, \qquad q = e^{-2x}

which is 1 - small and cancels nothing. q rather than expm1(2x) because it never overflows (q -> 0 is the correct limit and L(inf) = 1 falls out), and it also gives the derivative for free: csch^2(x) = 4q/(1-q)^2.

Inverse

L^-1 has a simple pole at y = 1, and near it x = 1/(1-y) - 2x^2 e^{-2x}, so for y >= 0.85 the seed is 1/(1-y) itself (relative error 2.2e-5 at 0.85, 4e-8 at 0.9, below u past 0.95). Below 0.85 the seed is y * q(y^2) / (1 - y^2) with q a minimax fit of L^-1(y)(1-y^2)/y, the same shape as Cohen’s Pade (3-y^2)/(1-y^2), which is what the vMF literature calls the Banerjee estimator. Both seeds share the one division 1/(1-y^2), since 1/(1-y) = (1+y)/(1-y^2).

The seed is then polished with one step. f32 takes Newton, x <- x - f/f' for f = L(x) - y: the error squares with constant ~1, so the 8e-5 seed lands at ~6e-9, past f32. f64 takes Halley, x <- x - 2ff'/(2f'^2 - ff''): L'' is nearly free on both branches (2 csch^2 coth - 2/x^3 from the same q, and the differentiated identity below 2), the step still has one division, and the error cubes with constant f'''/(6f') - (f''/(2f'))^2 at most ~0.07 (small y) and ~2e-4 at the tail crossover. So f64’s 1.1e-6 seed (deg 8 rather than f32’s deg 4 exactly for this) reaches full precision in one step where Newton needed two, i.e. a second exp and division. L is monotone and concave on x > 0, so from any positive seed either iteration is safe without safeguards.

The residual is formed as ((1-y) - 1/x) + 2q/(1-q) on the large branch rather than L(x) - y: 1 - y is exact for y >= 0.5, and that keeps the step accurate to u even where L(x) is within an ulp of 1, which L(x) - y cannot do (its error, u, divided by L' ~ 1/x^2, would grow as u x).

L^-1 itself is ill-conditioned near 1 (a relative error u in y moves the result by u/(1-y)), so callers with 1 - y in hand should compute it exactly before rounding, exactly as they would for acos near 1.

Functions

inv_langevin

L^-1(y) (or, with ONE_MINUS, L^-1(1 - t) from t directly): seed plus refine_steps of Newton (HALLEY = false) or Halley, see the module docs. Halley’s L'' costs a reciprocal and a few FMAs on top of Newton, which buys f64 a whole second step. f32 is already done after one Newton and would only pay.

langevin_primal

L(x) (or, with ONE_MINUS, 1 - L(x)) and L'(x) together.

refine_steps

Refinement steps for a given precision policy. One step (Newton from f32’s ~8e-5 seed, Halley from f64’s ~1e-6 seed, see the module docs) reaches the type’s full precision. Reference takes a second for good measure, Worst ships the seed.