lambert_w

This crate provides fast and accurate evaluation of the real valued parts of the principal and secondary branches of the Lambert W function with the method of Toshio Fukushima [1].

This method does not allocate, recurse, or iterate. It works by dividing the function's domain into subdomains. On each one, it uses a simple transformation of the input inserted into a rational function to approximate the true value.

The implementation uses conditional switches on the input value to select the appropriate subdomain, followed by either a square root (and possibly a division) or a logarithm. Then it performs a series of additions and multiplications by constants from a look-up table, and finishes the calculation with a division.

This crate provides two approximations of each branch, one with 50 bits of accuracy (implemented on 64-bit floats) and one with 24 bits (implemented on 32- and 64-bit floats). The one with 50 bits of accuracy uses higher degree polynomials in the rational functions compared to the one with only 24 bits, and thus more of the multiplications and additions by constants.

This crate can evaluate the approximation with 24 bits of accuracy on 32-bit floats, even though it is defined on 64-bit floats in Fukushima's paper. This may result in a reduction in the accuracy to less than 24 bits, but this reduction has not been quantified by the author of this crate.

This crate is no_std compatible, but can optionally depend on the standard library through features for a potential performance gain.

Examples

Compute the value of the omega constant with the principal branch of the Lambert W function:

use lambert_w::lambert_w0;
use approx::assert_abs_diff_eq;

let Ω = lambert_w0(1.0);

assert_abs_diff_eq!(Ω, 0.5671432904097839);

Evaluate the secondary branch of the Lambert W function at -ln(2)/2:

use lambert_w::lambert_wm1;
use approx::assert_abs_diff_eq;

let mln4 = lambert_wm1(-f64::ln(2.0) / 2.0);

assert_abs_diff_eq!(mln4, -f64::ln(4.0));

Do it on 32-bit floats:

use lambert_w::{lambert_w0f, lambert_wm1f};
use approx::assert_abs_diff_eq;

let Ω = lambert_w0f(1.0);
let mln4 = lambert_wm1f(-f32::ln(2.0) / 2.0);

assert_abs_diff_eq!(Ω, 0.56714329);
assert_abs_diff_eq!(mln4, -f32::ln(4.0));

The implementation can handle extreme inputs just as well:

use lambert_w::{lambert_w0, lambert_wm1};
use approx::assert_relative_eq;

let big = lambert_w0(f64::MAX);
let tiny = lambert_wm1(-f64::MIN_POSITIVE);

assert_relative_eq!(
    big,
    703.2270331047702,
    // Since the approximation used in this
    // example is accurate to 50 bits
    // it will sometimes have an
    // error larger than f64::EPSILON.
    max_relative = 1.5 * f64::EPSILON
);
assert_relative_eq!(tiny, -714.9686572379665);

Importing the LambertW trait lets you call the functions with postfix notation:

use lambert_w::LambertW;
use approx::assert_abs_diff_eq;

let ln1k = (1000.0 * f64::ln(1000.0)).lambert_w0();

assert_abs_diff_eq!(ln1k, f64::ln(1000.0));

The macros in the examples above are from the approx crate, and are used in the documentation examples of this crate. The assertion passes if the two supplied values are the same to within floating point error, or within an optional epsilon or relative difference.

Features

One of the below features must be enabled:

libm (enabled by default): if the std feature is disabled, this feature uses the libm crate to compute square roots and logarithms during function evaluation instead of the standard library.

std: use the standard library to compute square roots and logarithms for a potential performance gain. When this feature is disabled the crate is no_std compatible.

References

[1]: Toshio Fukushima. Precise and fast computation of Lambert W function by piecewise minimax rational function approximation with variable transformation. DOI: 10.13140/RG.2.2.30264.37128. November 2020.

⬆️ Back to top.

License