Langevin Function: Accurate Evaluation

Compute accurate values of the Langevin function (coth(x)-1/x) over its entire domain (-inf<x<inf).

The Langevin function arises frequently in the mathematical physics of electromagnetism. It is an entire function, defined by a simple formula:

L(x) = coth(x) - 1/x ,

with a limiting value of 0 at x=0.

The function is mathematically well conditioned, but the defining formula is not well suited for numerical evaluation for 'x' close to 0. A rapidly converging, and numerically stable, continued-fraction approximation is utlized instead for 0<|x|<1.

Numerical benchmarking suggests that the double-precision relative accuracy is roughly 15 decimal digits over all of -inf<x<inf.