# Accurate and efficient computation of the inverse Langevin function

The Langevin function $$\mathcal{L}(x) = \mathrm{coth}(x) - \frac{1}{x}$$ occurs in computations related to elastomers and paramagnetic materials. It is easily computed accurately and with high performance as long as the issue of subtractive cancellation near the origin is addressed, e.g. by using a minimax polynomial approximation in that region.

To eliminate the cancellation issue completely, choosing a good switch-over point to the polynomial approximation is important; some implementations that switch over at unity, which is suboptimal. A switchover point closer to $$2.0$$ is advisable. With a faithfully-rounded implementation of the exponential function $$\exp$$, which is usually provided by modern math libraries, a maximum error of $$1.3$$ ulp can be achieved in a straightforward manner. Exemplary C implementations are shown below.

However, many applications require the inverse Langevin function $$\mathcal{L}^{-1}(x)$$, which has no closed-form representation; numeric approximation is required. Most of the approximations used in the literature are quite inaccurate, providing the relative error in the range of 0.05% to 2%. Useful overviews are provided by:

Benjamin C. Marchi and Ellen M. Arruda, "An error-minimizing approach to inverse Langevin approximations." Rheol Acta (2015) 54:887-902

Martin Kröger, "Simple, admissible, and accurate approximants of the inverse Langevin and Brillouin functions, relevant for strong polymer deformations and flows." Journal of Non-Newtonian Fluid Mechanics 223 (2015) 77-87

Recent work shows that a lack of accuracy in the computation of the inverse Langevin function has a negative impact on the overall accuracy of the scientific computations in which it used:

Amine Ammar, "Effect of the inverse Langevin approximation on the solution of the Fokker-Planck equation of non-linear dilute polymer." Journal of Non-Newtonian Fluid Mechanics 231 (2016): 1-5.

The accurate computation of the inverse Langevin function is therefore highly desirable. In addition to accuracy, some use cases benefit from specific functional forms to ease of substitution, differentiation, and integration. Leaving aside such additional requirements, and focusing on minimizing relative error, how can the inverse Langevin function $$\mathcal{L}^{-1}(x)$$ be computed accurately (close to machine precision) and with good performance?

#include <math.h>  // fabs, exp, copysign

/* compute Langevin function L(a) = coth(a) - 1/a. Max error typically < 1.3 ulp */
float langevinf (float a)
{
float r, t;

t = fabsf (a);
if (t > 1.8125f) {
float e;
e = 1.0f / (expf (2.0f * t) - 1.0f);
r = (1.0f - 1.0f / t) + 2.0f * e;
r = copysignf (r, a);
} else {
float s;
s = a * a;
r =         7.70960469e-8f;
r = r * s - 1.65101926e-6f;
r = r * s + 2.03457112e-5f;
r = r * s - 2.10521728e-4f;
r = r * s + 2.11580913e-3f;
r = r * s - 2.22220998e-2f;
r = r * s + 8.33333284e-2f;
r = r * a + 0.25f * a;
}
return r;
}

/* compute Langevin function L(a) = coth(a) - 1/a. Max error typically < 1.3 ulp */
double langevin (double a)
{
double r, t;

t = fabs (a);
if (t > 1.8125) {
double e;
e = 1.0 / (exp (2.0 * t) - 1.0);
r = (1.0 - 1.0 / t) + 2.0 * e;
r = copysign (r, a);
} else {
double s, t, u;
s = a * a;
t = s * s;
r = 2.6224020345921745e-16 * s - 9.0382507125284313e-15;
u = 1.6058476918889783e-13 * s - 2.0544032928760129e-12;
r = r * t + u;
u = 2.2318295632954910e-11 * s - 2.2645716105162602e-10;
r = r * t + u;
u = 2.2484176131435534e-09 * s - 2.2212056151537571e-08;
r = r * t + u;
u = 2.1925751794205211e-07 * s - 2.1644032423757088e-06;
r = r * t + u;
u = 2.1377798795067705e-05 * s - 2.1164021156459166e-04;
r = r * t + u;
u = 2.1164021163937877e-03 * s - 2.2222222222221859e-02;
r = r * t + u;
r = r * s + 8.3333333333333329e-02;
r = r * a + 0.25 * a;
}
return r;
}


The inverse Langevin function $$\mathcal{L}^{-1}(x)$$ is an odd function. Therefore one needs to consider only approximation on the interval $$[0, 1]$$; the negative half-plane is treated via symmetry about the origin. Because of the singularity at unity, a polynomial approximation is not a fruitful direction of the investigation, although some authors have explored high-degree Taylor expansions around zero, for example

Mikhail Itskov, Roozbeh Dargazany, and Karl Hörnes, "Taylor expansion of the inverse function with application to the Langevin function." Mathematics and Mechanics of Solids 17, No. 7 (2012): 693-701

The 115-term series expansion from the paper is not practical due to the large amount of computation required. Rational approximations are well suited to approximating functions with singularities, and most work in the literature has focused on these, in particular in the form of Pade approximants. In recent years, the approach has been refined by adding compensation for the error terms, for example

Benjamin C. Marchi and Ellen M. Arruda. "Generalized error-minimizing, rational inverse Langevin approximations." Mathematics and Mechanics of Solids (2018): 1081286517754131.

Some authors have tried to approximate the shape of the inverse Langevin function with trigonometric functions, in particular, $$\tan$$ and $$\sin$$:

Jörgen Stefan Bergström, "Large strain time-dependent behavior of elastomeric materials." PhD diss., Massachusetts Institute of Technology, 1999.

Grant Keady, "The Langevin function and truncated exponential distributions." arXiv preprint arXiv:1501.02535 (2015).

Rafayel Petrosyan, "Improved approximations for some polymer extension models." Rheologica Acta 56, no. 1 (2017): 21-26.

While newer works push relative error down to about $$10^{-5}$$, I am aware of only one paper that tries to achieve close to machine precision. It does this by using piece-wise approximation with thousands of intervals to each of which it fits a cubic spline. The storage requirements for tables of coefficients reach to the hundreds of KBytes. This exceeds the size of first level caches on most processors, so the scheme may not be suitable in practice:

José María Benitez and Francisco Javier Montáns, "A simple and efficient numerical procedure to compute the inverse Langevin function with high accuracy". arXiv preprint arXiv:1806.08068 (2018).

For a different approach, I noted the similarity in overall shape between the inverse Langevin function and the inverse error function. As I pointed out in another answer, the latter can be approximated efficiently by making use of the transform $$t=\log(1-x^2)$$ and then generating approximations of the form $$x \cdot P(t)$$, where $$P$$ is a polynomial. Since $$\mathcal{L}(x) \approx 1 - \frac{1}{x}$$ near unity, it follows that $$\frac{1}{x-1}$$ is a good approximation to $$\mathcal{L}^{-1}(x)$$ near the singularity. Experimentally, for IEEE-754 single precision, using that approximation is suitable on $$[\frac{57}{64}, 1]$$, maintaining the maximum error of less than $$2.6$$ ulp.

My experiments show that, as expected, the approximation of $$\mathcal{L}^{-1}(x)$$ in the same manner as for the inverse error function requires just moderately more terms in the polynomial approximation. However, generating a high-quality minimax approximation requires some effort as the approximation problem is somewhat poorly conditioned. Starting with the Remez algorithms to generate an initial approximation, I used a heuristic refinement (similar to simulated annealing) to arrive at the final result. My preliminary experiments show that this approach is extensible for use cases requiring more than single-precision accuracy, but too computational expensive for reaching full double-precision accuracy, which would require polynomials of degree 40 or so.

Most modern hardware platforms provide a fused multiply-add operation (FMA) that helps reduce rounding errors and provides protection against some instances of subtractive cancellation. Its use also tends to improve performance, as is execution time is often identical or just slightly higher than a floating-point multiplication. The operation is often exposed through a library function, such as fma() in C and C++. However, this maps to slow emulation on hardware platforms without support for FMA, and some emulations have functional bugs.

An exemplary C implementation for IEEE-754 (2008) single precision (binary32) is shown below. While the accuracy achieved will depend on the quality of the standard math library's implementation of $$\log$$, most modern math libraries provide a faithfully-rounded implementation. With the three different implementations of logf, I tried the error bounds stated in the code comments were maintained.

#include <math.h> // fabs, log, copysign, fma

/*
Compute inverse Langevin function accurate to close to machine precision

USE_FMA == 0: max. ulp error < 4.37, max. relative error < 4.6e-7
USE_FMA == 1: max. ulp error < 3.63, max. relative error < 3.9e-7
*/
float langevininvf (float x)
{
float p, r, t;
if ((fabsf (x) > 0.890625f) && (fabsf (x) <= 1.0f)) {
r = copysignf (1.0f / (fabsf (x) - 1.0f), x);
} else {
#if USE_FMA
t = fmaf (x, 0.0f - x, 1.0f); // compute 1-x*x accurately
t = logf (t);
p =              2.18808651e-4f;  //  0x1.cae000p-13
p = fmaf (p, t, -7.90075958e-3f); // -0x1.02e460p-7
p = fmaf (p, t, -7.12909475e-2f); // -0x1.2401fap-4
p = fmaf (p, t, -2.40409255e-1f); // -0x1.ec5bb0p-3
p = fmaf (p, t, -4.14386481e-1f); // -0x1.a854eep-2
p = fmaf (p, t, -4.05751973e-1f); // -0x1.9f7d72p-2
p = fmaf (p, t, -2.56382823e-1f); // -0x1.068938p-2
p = fmaf (p, t, -1.22061908e-1f); // -0x1.f3f730p-4
p = fmaf (p, t,  5.00488281e-2f); //  0x1.9a0000p-5
p = fmaf (p, t, -1.84208527e-1f); // -0x1.794252p-3
p = fmaf (p, t,  3.98338169e-1f); //  0x1.97e5f6p-2
p = fmaf (p, t, -9.00006175e-1f); // -0x1.cccd9cp-1
p = fmaf (p, t,  4.99999970e-1f); //  0x1.fffffep-2
t = x + x;
r = fmaf (p, t, t);
#else // USE_FMA
t = (1.0f - x) * (1.0f + x); // compute 1-x*x accurately
t = logf (t);
p =         2.18808651e-4f; //  0x1.cae000p-13
p = p * t - 7.90080987e-3f; // -0x1.02e4ccp-7
p = p * t - 7.12910667e-2f; // -0x1.24021ap-4
p = p * t - 2.40408912e-1f; // -0x1.ec5b82p-3
p = p * t - 4.14385200e-1f; // -0x1.a85498p-2
p = p * t - 4.05750751e-1f; // -0x1.9f7d20p-2
p = p * t - 2.56382376e-1f; // -0x1.06891ap-2
p = p * t - 1.22061819e-1f; // -0x1.f3f718p-4
p = p * t + 5.00487499e-2f; //  0x1.99ffd6p-5
p = p * t - 1.84208632e-1f; // -0x1.794260p-3
p = p * t + 3.98337901e-1f; //  0x1.97e5e4p-2
p = p * t - 9.00006235e-1f; // -0x1.cccd9ep-1
p = p * t + 4.99999970e-1f; //  0x1.fffffep-2
t = x + x;
r = p * t + t;
#endif // USE_FMA
}
return r;
}


Computational platforms that provide fast hardware support for simple transcendental functions lend themselves to particularly efficient implementations of the proposed algorithm. The following shows an implementation for GPUs using CUDA that translates to only about twenty-five instructions.

/* max. ulp error < 3.93, max. relative error < 4.11e-7 */
__device__ float langevininvf (float x)
{
float fa, p, r, t;
fa = fabsf (x);
if ((fa > (57.0f / 64.0f)) && (fa <= 1.0f)) {
t = fa - 1.0f;
asm ("rcp.approx.ftz.f32 %0,%0;" : "+f"(t));
r = copysignf (t, x);
} else {
t = fmaf (x, -x, 1.0f);
asm ("lg2.approx.ftz.f32 %0,%0;" : "+f"(t));
p =              2.69152224e-6f;  //  0x1.694000p-19
p = fmaf (p, t, -1.40199525e-4f); // -0x1.26050cp-13
p = fmaf (p, t, -1.82510808e-3f); // -0x1.de70eep-10
p = fmaf (p, t, -8.87932442e-3f); // -0x1.22f52cp-7
p = fmaf (p, t, -2.20804978e-2f); // -0x1.69c452p-6
p = fmaf (p, t, -3.11916713e-2f); // -0x1.ff0b5ap-6
p = fmaf (p, t, -2.84342300e-2f); // -0x1.d1ddcep-6
p = fmaf (p, t, -1.95302274e-2f); // -0x1.3ffbb6p-6
p = fmaf (p, t,  1.15530221e-2f); //  0x1.7a91c6p-7
p = fmaf (p, t, -6.13459945e-2f); // -0x1.f68be0p-5
p = fmaf (p, t,  1.91382647e-1f); //  0x1.87f3a0p-3
p = fmaf (p, t, -6.23836875e-1f); // -0x1.3f678cp-1
p = fmaf (p, t,  5.00000000e-1f); //  0x1.000000p-1
t = x + x;
r = fmaf (p, t, t);
}
return r;
}