8
$\begingroup$

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;
}
Improve this question
$\endgroup$

1 Answer 1

Reset to default
12
$\begingroup$

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 showed 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. A rational approximation, as shown at the end of this answer, is more appropriate for a double-precision implementation.

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 almost machine precision

  USE_FMA == 0: max. ulp error < 4.27, max. relative error < 4.43e-7
  USE_FMA == 1: max. ulp error < 3.64, max. relative error < 3.84e-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.90076610e-3f); // -0x1.02e46ep-7
        p = fmaf (p, t, -7.12909698e-2f); // -0x1.240200p-4
        p = fmaf (p, t, -2.40409270e-1f); // -0x1.ec5bb2p-3
        p = fmaf (p, t, -4.14386481e-1f); // -0x1.a854eep-2
        p = fmaf (p, t, -4.05752033e-1f); // -0x1.9f7d76p-2
        p = fmaf (p, t, -2.56382942e-1f); // -0x1.068940p-2
        p = fmaf (p, t, -1.22061931e-1f); // -0x1.f3f736p-4
        p = fmaf (p, t,  5.00488468e-2f); //  0x1.9a000ap-5
        p = fmaf (p, t, -1.84208602e-1f); // -0x1.79425cp-3
        p = fmaf (p, t,  3.98338169e-1f); //  0x1.97e5f6p-2 
        p = fmaf (p, t, -9.00006115e-1f); // -0x1.cccd9ap-1 
        p = fmaf (p, t,  5.00000000e-1f); //  0x1.000000p-1
        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.18868256e-4f; //  0x1.cb0000p-13
        p = p * t - 7.90085457e-3f; // -0x1.02e52cp-7
        p = p * t - 7.12914243e-2f; // -0x1.24027ap-4
        p = p * t - 2.40408897e-1f; // -0x1.ec5b80p-3
        p = p * t - 4.14384902e-1f; // -0x1.a85484p-2
        p = p * t - 4.05750901e-1f; // -0x1.9f7d2ap-2
        p = p * t - 2.56382436e-1f; // -0x1.06891ep-2
        p = p * t - 1.22061789e-1f; // -0x1.f3f710p-4
        p = p * t + 5.00487871e-2f; //  0x1.99ffeap-5
        p = p * t - 1.84208602e-1f; // -0x1.79425cp-3
        p = p * t + 3.98338020e-1f; //  0x1.97e5ecp-2
        p = p * t - 9.00006115e-1f; // -0x1.cccd9ap-1
        p = p * t + 5.00000000e-1f; //  0x1.000000p-1
        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.925, max. relative error < 3.95e-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.40199758e-4f); // -0x1.26052cp-13
        p = fmaf (p, t, -1.82510854e-3f); // -0x1.de70f6p-10
        p = fmaf (p, t, -8.87932349e-3f); // -0x1.22f52ap-7
        p = fmaf (p, t, -2.20804960e-2f); // -0x1.69c450p-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;
}

As noted earlier, a rational approximation appears the best approach for a double-precision implementation. Noting the similarity of $\mathcal{L}^{-1}(x)$ to $\mathrm{tanh}(x)$, one approach is to transform the argument $x$ into $t=\frac{\ln(1+x)}{\ln(1-x)}$, then compute $\mathcal{L}^{-1}(x) \approx xR_{m,n}(t^2)$, where $R_{m,n}$ is a rational minimax approximation comprising a polynomial $p$ of degree $m$ in the numerator and a polynomial $q$ of degree $n$ in the denominator. Trying to numerically compute the coefficients of the minimax approximation with a variant of the Remez algorithm turns out to be extremely ill conditioned. When increasing $m, n$ to achieve full double-precision accuracy, the solutions tend to degenerate: the number of peaks in the error curve is deficient and the equioscillation property is lost. A recent paper reports similar difficulties in approximating $\mathcal{L}^{-1}(x)$ on $[0,1)$:

Nicolas Brisebarre and Silviu-Ioan Filip, "Towards Machine-Efficient Rational $L^∞$-Approximations of Mathematical Functions," Draft publication hal-04093020, May 2023 (online).

The best I have managed so far is to create an approximation $R_{14,14}$ that results in $\mathcal{L}^{-1}(x)$ being computed with a maximum relative error of $6 \cdot 10^{-15}$ or, using a different error metric, 47 double-precision ulp. These error bounds are based on 200M random test vectors. My attempts of creating an approximation $R_{15,15}$ have failed so far, with minute changes in the starting conditions of the Remez algorithm leading to vastly different sets of polynomial coefficients, none of them representing a non-degenerate solution.

/* Compute inverse Langevin function with maximum relative error of 6e-15 */
double langevininv (double x)
{
    double fa, p, q, r, t;
    fa = fabs (x);
    if ((fa > 971.0/1024.0) && (fa <= 1.0)) {
        // maximum error: 4 ulps
        t = fa - 1.0;
        t = 1.0 / t;
        r = copysign (t, x);
    } else {
        // maximum error: 47 ulps 
        t = log1p (2.0 * fa / (1.0 - fa));
        t = t * t;
        p =             8.7442384822583792e-14;  //  0x1.89ce34b44b32ep-44
        p = fma (p, t, -1.0246877312258876e-12); // -0x1.206c8855594fap-40
        p = fma (p, t,  6.9112635465060234e-11); //  0x1.2ff5e8eabbf42p-34
        p = fma (p, t, -1.5954120724949512e-9);  // -0x1.b68b24dc5af04p-30
        p = fma (p, t,  3.7719321864398964e-8);  //  0x1.4401aa7c94e3bp-25
        p = fma (p, t, -5.9501169268478895e-7);  // -0x1.3f71c8cae9a0ap-21
        p = fma (p, t,  7.8764836653194943e-6);  //  0x1.084a7ac083ecfp-17
        p = fma (p, t, -7.3715838816829616e-5);  // -0x1.352fc77b8b1c1p-14
        p = fma (p, t,  5.3389955343705986e-4);  //  0x1.17eac8f88c356p-11
        p = fma (p, t, -2.3748559515852675e-3);  // -0x1.3746f14a1f408p-9 
        p = fma (p, t,  5.7661901230895048e-3);  //  0x1.79e49e004dd54p-8
        p = fma (p, t,  1.6354150376575827e-2);  //  0x1.0bf247415c4fdp-6
        p = fma (p, t, -1.0764052701852338e-1);  // -0x1.b8e545f44cd9bp-4
        p = fma (p, t,  2.9732462670549570e-1);  //  0x1.3075ddeffb679p-2
        p = fma (p, t,  1.0000000000000000e+0);  //  0x1.0000000000000p+0
        q =            -1.8576261900157100e-16;  // -0x1.ac56e54dc47acp-53
        q = fma (q, t,  3.7025672616258288e-14); //  0x1.4d7f636bcc935p-45
        q = fma (q, t, -3.8488081523689396e-12); // -0x1.0ed5f6be6f859p-38
        q = fma (q, t,  2.3670115378123239e-10); //  0x1.044173a5c110cp-32
        q = fma (q, t, -7.2533635216156031e-9);  // -0x1.f2728540cd2dcp-28
        q = fma (q, t,  1.6722199004755106e-7);  //  0x1.671b42e087ca4p-23
        q = fma (q, t, -2.8147043983109852e-6);  // -0x1.79c881b7b21b0p-19
        q = fma (q, t,  3.6736019004458165e-5);  //  0x1.342a0006dc6d1p-15
        q = fma (q, t, -3.6078759090487885e-4);  // -0x1.7a502e75b71c8p-12
        q = fma (q, t,  2.6901070708282545e-3);  //  0x1.609903c8ad1bap-9
        q = fma (q, t, -1.4214262366936270e-2);  // -0x1.d1c5e0005a453p-7
        q = fma (q, t,  4.9851609255380437e-2);  //  0x1.986266ecf45ddp-5
        q = fma (q, t, -7.0008037607361850e-2);  // -0x1.1ec0bf7fb12e8p-4
        q = fma (q, t, -1.5267537329450756e-1);  // -0x1.38adddb9a2d1dp-3
        q = fma (q, t,  1.0000000000000000e+0);  //  0x1.0000000000000p+0
        r = fma (p / q, x, x + x);
    }
    return r;
}
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.