Boost C++ Libraries

Synopsis

#include <boost/math/special_functions/lambert_w.hpp>

namespace boost { namespace math {

  template <class T>
  calculated-result-type lambert_w0(T z);                        // W0 branch, default policy.
  template <class T>
  calculated-result-type lambert_wm1(T z);                       // W-1 branch, default policy.
  template <class T>
  calculated-result-type lambert_w0_prime(T z);                  // W0 branch 1st derivative.
  template <class T>
  calculated-result-type lambert_wm1_prime(T z);                 // W-1 branch 1st derivative.

  template <class T, class Policy>
  calculated-result-type lambert_w0(T z, const Policy&);         // W0 with policy.
  template <class T, class Policy>
  calculated-result-type lambert_wm1(T z, const Policy&);        // W-1 with policy.
  template <class T, class Policy>
  calculated-result-type lambert_w0_prime(T z, const Policy&);   // W0 derivative with policy.
  template <class T, class Policy>
  calculated-result-type lambert_wm1_prime(T z, const Policy&);  // W-1 derivative with policy.

 } // namespace boost
 } // namespace math

Description

The Lambert W function is the solution of the equation W(z)e^W(z) = z. It is also called the Omega function, the inverse of f(W) = We^W.

On the interval [0, ∞), there is just one real solution. On the interval (-e^-1, 0), there are two real solutions, generating two branches which we will denote by W₀ and W_-1. In Boost.Math, we call these principal branches lambert_w0 and lambert_wm1; their derivatives are labelled lambert_w0_prime and lambert_wm1_prime.

There is a singularity where the branches meet at e^-1 ≅ -0.367879. Approaching this point, the condition number of function evaluation tends to infinity, and the only method of recovering high accuracy is use of higher precision.

This implementation computes the two real branches W₀ and W_-1 with the functions lambert_w0 and lambert_wm1, and their derivatives, lambert_w0_prime and lambert_wm1_prime. Complex arguments are not supported.

The final Policy argument is optional and can be used to control how the function deals with errors. Refer to Policies for more details and see examples below.

Applications of the Lambert W function

The Lambert W function has a myriad of applications. Corless et al. provide a summary of applications, from the mathematical, like iterated exponentiation and asymptotic roots of trinomials, to the real-world, such as the range of a jet plane, enzyme kinetics, water movement in soil, epidemics, and diode current (an example replicated here). Since the publication of their landmark paper, there have been many more applications, and also many new implementations of the function, upon which this implementation builds.

Examples

The most basic usage of the Lambert-W function is demonstrated below:

#include <boost/math/special_functions/lambert_w.hpp> // For lambert_w function.

using boost::math::lambert_w0;
using boost::math::lambert_wm1;

std::cout.precision(std::numeric_limits<double>::max_digits10);
// Show all potentially significant decimal digits,
std::cout << std::showpoint << std::endl;
// and show significant trailing zeros too.

double z = 10.;
double r = lambert_w0(z); // Default policy for double.
std::cout << "lambert_w0(z) = " << r << std::endl;
// lambert_w0(z) = 1.7455280027406994

Other floating-point types can be used too, here float, including user-defined types like Boost.Multiprecision. It is convenient to use a function like show_value to display all (and only) potentially significant decimal digits, including any significant trailing zeros, (std::numeric_limits<T>::max_digits10) for the type T.

float z = 10.F;
float r;
r = lambert_w0(z);        // Default policy digits10 = 7, digits2 = 24
std::cout << "lambert_w0(";
show_value(z);
std::cout << ") = ";
show_value(r);
std::cout << std::endl;   // lambert_w0(10.0000000) = 1.74552798

Example of an integer argument to lambert_w0, showing that an int literal is correctly promoted to a double.

std::cout.precision(std::numeric_limits<double>::max_digits10);
double r = lambert_w0(10);                           // Pass an int argument "10" that should be promoted to double argument.
std::cout << "lambert_w0(10) = " << r << std::endl;  // lambert_w0(10) = 1.7455280027406994
double rp = lambert_w0(10);
std::cout << "lambert_w0(10) = " << rp << std::endl;
// lambert_w0(10) = 1.7455280027406994
auto rr = lambert_w0(10);                            // C++11 needed.
std::cout << "lambert_w0(10) = " << rr << std::endl;
// lambert_w0(10) = 1.7455280027406994 too, showing that rr has been promoted to double.

Using Boost.Multiprecision types to get much higher precision is painless.

cpp_dec_float_50 z("10");
// Note construction using a decimal digit string "10", 
// NOT a floating-point double literal 10.
cpp_dec_float_50 r;
r = lambert_w0(z);
std::cout << "lambert_w0("; show_value(z); std::cout << ") = ";
show_value(r);
std::cout << std::endl;
// lambert_w0(10.000000000000000000000000000000000000000000000000000000000000000000000000000000) = 
//   1.7455280027406993830743012648753899115352881290809413313533156980404446940000000

	Warning
	When using multiprecision, take very great care not to construct or assign non-integers from double, float ... silently losing precision. Use "1.2345678901234567890123456789" rather than 1.2345678901234567890123456789.

Using multiprecision types, it is all too easy to get multiprecision precision wrong!

cpp_dec_float_50 z(0.7777777777777777777777777777777777777777777777777777777777777777777777777);
// Compiler evaluates the nearest double-precision binary representation,
// from the max_digits10 of the floating_point literal double 0.7777777777777777777777777777...,
// so any extra digits in the multiprecision type 
// beyond max_digits10 (usually 17) are random and meaningless.
cpp_dec_float_50 r;
r = lambert_w0(z);
std::cout << "lambert_w0(";
show_value(z);
std::cout << ") = "; show_value(r);
std::cout << std::endl;
// lambert_w0(0.77777777777777779011358916250173933804035186767578125000000000000000000000000000)
//   = 0.48086152073210493501934682309060873341910109230469724725005039758139532631901386

	Note
	See spurious non-seven decimal digits appearing after digit #17 in the argument 0.7777777777777777...!

And similarly constructing from a literal double 0.9, with more random digits after digit number 17.

cpp_dec_float_50 z(0.9); // Construct from floating_point literal double 0.9.
cpp_dec_float_50 r;
r = lambert_w0(0.9);
std::cout << "lambert_w0(";
show_value(z);
std::cout << ") = "; show_value(r);
std::cout << std::endl;
// lambert_w0(0.90000000000000002220446049250313080847263336181640625000000000000000000000000000)
//   = 0.52983296563343440510607251781038939952850341796875000000000000000000000000000000
std::cout << "lambert_w0(0.9) = " << lambert_w0(static_cast<double>(0.9))
// lambert_w0(0.9) 
//   = 0.52983296563343441
  << std::endl;

Note how the cpp_float_dec_50 result is only as correct as from a double = 0.9.

Now see the correct result for all 50 decimal digits constructing from a decimal digit string "0.9":

cpp_dec_float_50 z("0.9");     // Construct from decimal digit string.
cpp_dec_float_50 r;
r = lambert_w0(z);
std::cout << "lambert_w0(";
show_value(z);
std::cout << ") = "; show_value(r);
std::cout << std::endl;
// 0.90000000000000000000000000000000000000000000000000000000000000000000000000000000)
// = 0.52983296563343441213336643954546304857788132269804249284012528304239956413801252

Note the expected zeros for all places up to 50 - and the correct Lambert W result!

(It is just as easy to compute even higher precisions, at least to thousands of decimal digits, but not shown here for brevity. See lambert_w_simple_examples.cpp for comparison of an evaluation at 1000 decimal digit precision with Wolfram Alpha).

Policies can be used to control what action to take on errors:

// Define an error handling policy:
typedef policy<
  domain_error<throw_on_error>,
  overflow_error<ignore_error> // possibly unwise?
> my_throw_policy;

std::cout.precision(std::numeric_limits<double>::max_digits10);
// Show all potentially significant decimal digits,
std::cout << std::showpoint << std::endl;
// and show significant trailing zeros too.
double z = +1;
std::cout << "Lambert W (" << z << ") = " << lambert_w0(z) << std::endl;
// Lambert W (1.0000000000000000) = 0.56714329040978384
std::cout << "\nLambert W (" << z << ", my_throw_policy()) = "
  << lambert_w0(z, my_throw_policy()) << std::endl;
// Lambert W (1.0000000000000000, my_throw_policy()) = 0.56714329040978384

An example error message:

Error in function boost::math::lambert_wm1<RealType>(<RealType>):
Argument z = 1 is out of range (z <= 0) for Lambert W-1 branch! (Try Lambert W0 branch?)

Showing an error reported if a value is passed to lambert_w0 that is out of range, (and was probably meant to be passed to lambert_wm1 instead).

double z = +1.;
double r = lambert_wm1(z);
std::cout << "lambert_wm1(+1.) = " << r << std::endl;

The full source of these examples is at lambert_w_simple_examples.cpp

Diode Resistance Example

A typical example of a practical application is estimating the current flow through a diode with series resistance from a paper by Banwell and Jayakumar.

Having the Lambert W function available makes it simple to reproduce the plot in their paper (Fig 2) comparing estimates using with Lambert W function and some actual measurements. The colored curves show the effect of various series resistance on the current compared to an extrapolated line in grey with no internal (or external) resistance.

Two formulae relating the diode current and effect of series resistance can be combined, but yield an otherwise intractable equation relating the current versus voltage with a varying series resistance. This was reformulated as a generalized equation in terms of the Lambert W function:

Banwell and Jakaumar equation 5

    I(V) = μ V_T/ R _S ․ W₀(I₀ R_S / (μ V_T))

Using these variables

double nu = 1.0; // Assumed ideal.
double vt = v_thermal(25); // v thermal, Shockley equation, expect about 25 mV at room temperature.
double boltzmann_k = 1.38e-23; // joules/kelvin
double temp = 273 + 25;
double charge_q = 1.6e-19; // column
vt = boltzmann_k * temp / charge_q;
std::cout << "V thermal " << vt << std::endl; // V thermal 0.0257025 = 25 mV
double rsat = 0.;
double isat = 25.e-15; //  25 fA;
std::cout << "Isat = " << isat << std::endl;
double re = 0.3;  // Estimated from slope of straight section of graph (equation 6).
double v = 0.9;
double icalc = iv(v, vt, 249., re, isat);
std::cout << "voltage = " << v << ", current = " << icalc << ", " << log(icalc) << std::endl; // voltage = 0.9, current = 0.00108485, -6.82631

the formulas can be rendered in C++

double iv(double v, double vt, double rsat, double re, double isat, double nu = 1.)
{
  // V thermal 0.0257025 = 25 mV
  // was double i = (nu * vt/r) * lambert_w((i0 * r) / (nu * vt)); equ 5.

  rsat = rsat + re;
  double i = nu * vt / rsat;
 // std::cout << "nu * vt / rsat = " << i << std::endl; // 0.000103223

  double x = isat * rsat / (nu * vt);
//  std::cout << "isat * rsat / (nu * vt) = " << x << std::endl;

  double eterm = (v + isat * rsat) / (nu * vt);
 // std::cout << "(v + isat * rsat) / (nu * vt) = " << eterm << std::endl;

  double e = exp(eterm);
//  std::cout << "exp(eterm) = " << e << std::endl;

  double w0 = lambert_w0(x * e);
//  std::cout << "w0 = " << w0 << std::endl;
  return i * w0 - isat;
} // double iv

to reproduce their Fig 2:

The plotted points for no external series resistance (derived from their published plot as the raw data are not publicly available) are used to extrapolate back to estimate the intrinsic emitter resistance as 0.3 ohm. The effect of external series resistance is visible when the colored lines start to curve away from the straight line as voltage increases.

See lambert_w_diode.cpp and lambert_w_diode_graph.cpp for details of the calculation.

Existing implementations

The principal value of the Lambert W function is implemented in the Wolfram Language as ProductLog[k, z], where k is the branch.

The symbolic algebra program Maple also computes Lambert W to an arbitrary precision.

Controlling the compromise between Precision and Speed

Floating-point types double and float

This implementation provides good precision and excellent speed for __fundamental float and double.

All the functions usually return values within a few Unit in the last place (ULP) for the floating-point type, except for very small arguments very near zero, and for arguments very close to the singularity at the branch point.

By default, this implementation provides the best possible speed. Very slightly average higher precision and less bias might be obtained by adding a Halley step refinement, but at the cost of more than doubling the runtime.

Floating-point types larger than double

For floating-point types with precision greater than double and float fundamental (built-in) types, a double evaluation is used as a first approximation followed by Halley refinement, using a single step where it can be predicted that this will be sufficient, and only using Halley iteration when necessary. Higher precision types are always going to be very, very much slower.

The 'best' evaluation (the nearest representable) can be achieved by static_casting from a higher precision type, typically a Boost.Multiprecision type like cpp_bin_float_50, but at the cost of increasing run-time 100-fold; this has been used here to provide some of our reference values for testing.

For example, we get a reference value using a high precision type, for example;

using boost::multiprecision::cpp_bin_float_50;

that uses Halley iteration to refine until it is as precise as possible for this cpp_bin_float_50 type.

As a further check we can compare this with a Wolfram Alpha computation using command N[ProductLog[10.], 50] to get 50 decimal digits and similarly N[ProductLog[10.], 17] to get the nearest representable for 64-bit double precision.

 using boost::multiprecision::cpp_bin_float_50;
 using boost::math::float_distance;

 cpp_bin_float_50 z("10."); // Note use a decimal digit string, not a double 10.
 cpp_bin_float_50 r;
 std::cout.precision(std::numeric_limits<cpp_bin_float_50>::digits10);

 r = lambert_w0(z); // Default policy.
 std::cout << "lambert_w0(z) cpp_bin_float_50  = " << r << std::endl;
 //lambert_w0(z) cpp_bin_float_50  = 1.7455280027406993830743012648753899115352881290809
 //       [N[productlog[10], 50]] == 1.7455280027406993830743012648753899115352881290809
 std::cout.precision(std::numeric_limits<double>::max_digits10);
 std::cout << "lambert_w0(z) static_cast from cpp_bin_float_50  = "
   << static_cast<double>(r) << std::endl;
 // double lambert_w0(z) static_cast from cpp_bin_float_50  = 1.7455280027406994
 // [N[productlog[10], 17]]                                == 1.7455280027406994
std::cout << "bits different from Wolfram = "
  << static_cast<int>(float_distance(static_cast<double>(r), 1.7455280027406994))
  << std::endl; // 0

giving us the same nearest representable using 64-bit double as 1.7455280027406994.

However, the rational polynomial and Fukushima Schroder approximations are so good for type float and double that negligible improvement is gained from a double Halley step.

This is shown with lambert_w_precision_example.cpp for Lambert W₀:

using boost::math::lambert_w_detail::lambert_w_halley_step;
using boost::math::epsilon_difference;
using boost::math::relative_difference;

std::cout << std::showpoint << std::endl; // and show any significant trailing zeros too.
std::cout.precision(std::numeric_limits<double>::max_digits10); // 17 decimal digits for double.

cpp_bin_float_50 z50("1.23"); // Note: use a decimal digit string, not a double 1.23!
double z = static_cast<double>(z50);
cpp_bin_float_50 w50;
w50 = lambert_w0(z50);
std::cout.precision(std::numeric_limits<cpp_bin_float_50>::max_digits10); // 50 decimal digits.
std::cout << "Reference Lambert W (" << z << ") =\n                                              "
  << w50 << std::endl;
std::cout.precision(std::numeric_limits<double>::max_digits10); // 17 decimal digits for double.
double wr = static_cast<double>(w50);
std::cout << "Reference Lambert W (" << z << ") =    " << wr << std::endl;

double w = lambert_w0(z);
std::cout << "Rat/poly Lambert W  (" << z << ")  =   " << lambert_w0(z) << std::endl;
// Add a Halley step to the value obtained from rational polynomial approximation.
double ww = lambert_w_halley_step(lambert_w0(z), z);
std::cout << "Halley Step Lambert W (" << z << ") =  " << lambert_w_halley_step(lambert_w0(z), z) << std::endl;

std::cout << "absolute difference from Halley step = " << w - ww << std::endl;
std::cout << "relative difference from Halley step = " << relative_difference(w, ww) << std::endl;
std::cout << "epsilon difference from Halley step  = " << epsilon_difference(w, ww) << std::endl;
std::cout << "epsilon for float =                    " << std::numeric_limits<double>::epsilon() << std::endl;
std::cout << "bits different from Halley step  =     " << static_cast<int>(float_distance(w, ww)) << std::endl;

with this output:

Reference Lambert W (1.2299999999999999822364316059974953532218933105468750) =
0.64520356959320237759035605255334853830173300262666480
Reference Lambert W (1.2300000000000000) =    0.64520356959320235
Rat/poly Lambert W  (1.2300000000000000)  =   0.64520356959320224
Halley Step Lambert W (1.2300000000000000) =  0.64520356959320235
absolute difference from Halley step = -1.1102230246251565e-16
relative difference from Halley step = 1.7207329236029286e-16
epsilon difference from Halley step  = 0.77494921535422934
epsilon for float =                    2.2204460492503131e-16
bits different from Halley step  =     1

and then for W_-1:

using boost::math::lambert_w_detail::lambert_w_halley_step;
using boost::math::epsilon_difference;
using boost::math::relative_difference;

std::cout << std::showpoint << std::endl; // and show any significant trailing zeros too.
std::cout.precision(std::numeric_limits<double>::max_digits10); // 17 decimal digits for double.

cpp_bin_float_50 z50("-0.123"); // Note: use a decimal digit string, not a double -1.234!
double z = static_cast<double>(z50);
cpp_bin_float_50 wm1_50;
wm1_50 = lambert_wm1(z50);
std::cout.precision(std::numeric_limits<cpp_bin_float_50>::max_digits10); // 50 decimal digits.
std::cout << "Reference Lambert W-1 (" << z << ") =\n                                                  "
  << wm1_50 << std::endl;
std::cout.precision(std::numeric_limits<double>::max_digits10); // 17 decimal digits for double.
double wr = static_cast<double>(wm1_50);
std::cout << "Reference Lambert W-1 (" << z << ") =    " << wr << std::endl;

double w = lambert_wm1(z);
std::cout << "Rat/poly Lambert W-1 (" << z << ")  =    " << lambert_wm1(z) << std::endl;
// Add a Halley step to the value obtained from rational polynomial approximation.
double ww = lambert_w_halley_step(lambert_wm1(z), z);
std::cout << "Halley Step Lambert W (" << z << ") =    " << lambert_w_halley_step(lambert_wm1(z), z) << std::endl;

std::cout << "absolute difference from Halley step = " << w - ww << std::endl;
std::cout << "relative difference from Halley step = " << relative_difference(w, ww) << std::endl;
std::cout << "epsilon difference from Halley step  = " << epsilon_difference(w, ww) << std::endl;
std::cout << "epsilon for float =                    " << std::numeric_limits<double>::epsilon() << std::endl;
std::cout << "bits different from Halley step  =     " << static_cast<int>(float_distance(w, ww)) << std::endl;

with this output:

Reference Lambert W-1 (-0.12299999999999999822364316059974953532218933105468750) =
-3.2849102557740360179084675531714935199110302996513384
Reference Lambert W-1 (-0.12300000000000000) =    -3.2849102557740362
Rat/poly Lambert W-1 (-0.12300000000000000)  =    -3.2849102557740357
Halley Step Lambert W (-0.12300000000000000) =    -3.2849102557740362
absolute difference from Halley step = 4.4408920985006262e-16
relative difference from Halley step = 1.3519066740696092e-16
epsilon difference from Halley step  = 0.60884463935795785
epsilon for float =                    2.2204460492503131e-16
bits different from Halley step  =     -1

Distribution of differences from 'best' double evaluations

The distribution of differences from 'best' are shown in these graphs comparing double precision evaluations with reference 'best' z50 evaluations using cpp_bin_float_50 type reduced to double with static_cast<double(z50) :

As noted in the implementation section, the distribution of these differences is somewhat biased for Lambert W_-1 and this might be reduced using a double Halley step at small runtime cost. But if you are seriously concerned to get really precise computations, the only way is using a higher precision type and then reduce to the desired type. Fortunately, Boost.Multiprecision makes this very easy to program, if much slower.

Edge and Corner cases

The W₀ Branch

The domain of W₀ is [-e^-1, ∞). Numerically,

(An infinite argument probably indicates that something has already gone wrong, but if it is desired to return infinity, this case should be handled before calling lambert_w0).

W_-1 Branch

The domain of W_-1 is [-e^-1, 0). Numerically,

Denormalized values are not supported for Lambert W_-1 (because not all floating-point types denormalize), and anyway it only covers a tiny fraction of the range of possible z arguments values.

Compilers

The lambert_w.hpp code has been shown to work on most C++98 compilers. (Apart from requiring C++11 extensions for using of std::numeric_limits<>::max_digits10 in some diagnostics. Many old pre-c++11 compilers provide this extension but may require enabling to use, for example using b2/bjam the lambert_w examples use this command:

[ run lambert_w_basic_example.cpp  : : : [ requires cxx11_numeric_limits ] ]

See jamfile.v2.)

For details of which compilers are expected to work see lambert_w tests and examples in:
Boost Test Summary report for master branch (used for latest release)
Boost Test Summary report for latest developer branch.

As expected, debug mode is very much slower than release.

Diagnostics Macros

Several macros are provided to output diagnostic information (potentially much output). These can be statements, for example:

#define BOOST_MATH_INSTRUMENT_LAMBERT_W_TERMS

placed before the lambert_w include statement

#include <boost/math/special_functions/lambert_w.hpp>,

or defined on the project compile command-line: /DBOOST_MATH_INSTRUMENT_LAMBERT_W_TERMS,

or defined in a jamfile.v2: <define>BOOST_MATH_INSTRUMENT_LAMBERT_W_TERMS

// #define-able macros
BOOST_MATH_INSTRUMENT_LAMBERT_W_HALLEY                     // Halley refinement diagnostics.
BOOST_MATH_INSTRUMENT_LAMBERT_W_PRECISION                  // Precision.
BOOST_MATH_INSTRUMENT_LAMBERT_WM1                          // W1 branch diagnostics.
BOOST_MATH_INSTRUMENT_LAMBERT_WM1_HALLEY                   // Halley refinement diagnostics only for W-1 branch.
BOOST_MATH_INSTRUMENT_LAMBERT_WM1_TINY                     // K > 64, z > -1.0264389699511303e-26
BOOST_MATH_INSTRUMENT_LAMBERT_WM1_LOOKUP                   // Show results from W-1 lookup table.
BOOST_MATH_INSTRUMENT_LAMBERT_W_SCHROEDER                  // Schroeder refinement diagnostics.
BOOST_MATH_INSTRUMENT_LAMBERT_W_TERMS                      // Number of terms used for near-singularity series.
BOOST_MATH_INSTRUMENT_LAMBERT_W_SINGULARITY_SERIES         // Show evaluation of series near branch singularity.
BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES
BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES_ITERATIONS  // Show evaluation of series for small z.

Implementation

There are many previous implementations, each with increasing accuracy and/or speed. See references below.

For most of the range of z arguments, some initial approximation followed by a single refinement, often using Halley or similar method, gives a useful precision. For speed, several implementations avoid evaluation of a iteration test using the exponential function, estimating that a single refinement step will suffice, but these rarely get to the best result possible. To get a better precision, additional refinements, probably iterative, are needed for example, using Halley or Schröder methods.

For C++, the most precise results possible, closest to the nearest representable for the C++ type being used, it is usually necessary to use a higher precision type for intermediate computation, finally static-casting back to the smaller desired result type. This strategy is used by Maple and WolframAlpha, for example, using arbitrary precision arithmetic, and some of their high-precision values are used for testing this library. This method is also used to provide some Boost.Test values using Boost.Multiprecision, typically, a 50 decimal digit type like cpp_bin_float_50 static_cast to a float, double or long double type.

For z argument values near the singularity and near zero, other approximations may be used, possibly followed by refinement or increasing number of series terms until a desired precision is achieved. At extreme arguments near to zero or the singularity at the branch point, even this fails and the only method to achieve a really close result is to cast from a higher precision type.

In practical applications, the increased computation required (often towards a thousand-fold slower and requiring much additional code for Boost.Multiprecision) is not justified and the algorithms here do not implement this. But because the Boost.Lambert_W algorithms has been tested using Boost.Multiprecision, users who require this can always easily achieve the nearest representation for __fundamental - if the application justifies the very large extra computation cost.

Evolution of this implementation

One compact real-only implementation was based on an algorithm by Thomas Luu, Thesis, University College London (2015), (see routine 11 on page 98 for his Lambert W algorithm) and his Halley refinement is used iteratively when required. A first implementation was based on Thomas Luu's code posted at Boost Trac #11027. It has been implemented from Luu's algorithm but templated on RealType parameter and result and handles both fundamental (built-in) types (float, double, long double), Boost.Multiprecision, and also has been tested successfully with a proposed fixed_point type.

A first approximation was computed using the method of Barry et al (see references 5 & 6 below). This was extended to the widely used TOMS443 FORTRAN and C++ versions by John Burkardt using Schroeder refinement(s). (For users only requiring an accuracy of relative accuracy of 0.02%, Barry's function alone might suffice, but a better rational function approximation method has since been developed for this implementation).

We also considered using Newton-Raphson iteration method.

f(w) = w e^w -z = 0 // Luu equation 6.37
f'(w) = e^w (1 + w), Wolfram alpha (d)/(dw)(f(w) = w exp(w) - z) = e^w (w + 1)
if (f(w) / f'(w) -1 < tolerance
w1 = w0 - (expw0 * (w0 + 1)); // Refine new Newton/Raphson estimate.

but concluded that since the Newton-Raphson method takes typically 6 iterations to converge within tolerance, whereas Halley usually takes only 1 to 3 iterations to achieve an result within 1 Unit in the last place (ULP), so the Newton-Raphson method is unlikely to be quicker than the additional cost of computing the 2nd derivative for Halley's method.

Halley refinement uses the simplified formulae obtained from Wolfram Alpha

[2(z exp(z)-w) d/dx (z exp(z)-w)] / [2 (d/dx (z exp(z)-w))^2 - (z exp(z)-w) d^2/dx^2 (z exp(z)-w)]

Implementing Compact Algorithms

The most compact algorithm can probably be implemented using the log approximation of Corless et al. followed by Halley iteration (but is also slowest and least precise near zero and near the branch singularity).

Implementing Faster Algorithms

More recently, the Tosio Fukushima has developed an even faster algorithm, avoiding any transcendental function calls as these are necessarily expensive. The current implementation of Lambert W_-1 is based on his algorithm starting with a translation from Fukushima's FORTRAN into C++ by Darko Veberic.

Many applications of the Lambert W function make many repeated evaluations for Monte Carlo methods; for these applications speed is very important. Luu, and Chapeau-Blondeau and Monir provide typical usage examples.

Fukushima improves the important observation that much of the execution time of all previous iterative algorithms was spent evaluating transcendental functions, usually exp. He has put a lot of work into avoiding any slow transcendental functions by using lookup tables and bisection, finishing with a single Schroeder refinement, without any check on the final precision of the result (necessarily evaluating an expensive exponential).

Theoretical and practical tests confirm that Fukushima's algorithm gives Lambert W estimates with a known small error bound (several Unit in the last place (ULP)) over nearly all the range of z argument.

A mean difference was computed to express the typical error and is often about 0.5 epsilon, the theoretical minimum. Using the Boost.Math float_distance, we can also express this as the number of bits that are different from the nearest representable or 'exact' or 'best' value. The number and distribution of these few bits differences was studied by binning, including their sign. Bins for (signed) 0, 1, 2 and 3 and 4 bits proved suitable.

However, though these give results within a few machine epsilon of the nearest representable result, they do not get as close as is very often possible with further refinement, nrealy always to within one or two machine epsilon.

More significantly, the evaluations of the sum of all signed differences using the Fukshima algorithm show a slight bias, being more likely to be a bit or few below the nearest representation than above; bias might have unwanted effects on some statistical computations.

Fukushima's method also does not cover the full range of z arguments of 'float' precision and above.

For this implementation of Lambert W₀, John Maddock used the Boost.Math Remez algorithm method program to devise a rational function for several ranges of argument for the W₀ branch of Lambert W function. These minimax rational approximations are combined for an algorithm that is both smaller and faster.

Sadly it has not proved practical to use the same Remez algorithm method for Lambert W_-1 branch and so the Fukushima algorithm is retained for this branch.

An advantage of both minimax rational Remez algorithm approximations is that the distribution from the reference values is reasonably random and insignificantly biased.

For example, table below a test of Lambert W₀ 10000 values of argument covering the main range of possible values, 10000 comparisons from z = 0.0501 to 703, in 0.001 step factor 1.05 when module 7 == 0

Lambert W₀ values computed using the Fukushima method with Schroeder refinement gave about 1/6 lambert_w0 values that are one bit different from the 'best', and < 1% that are a few bits 'wrong'. If a Halley refinement step is added, only 1 in 30 are even one bit different, and only 2 two-bits 'wrong'.

With the rational polynomial approximation method, there are a third one-bit from the best and none more than two-bits. Adding a Halley step (or iteration) reduces the number that are one-bit different from about a third down to one in 30; this is unavoidable 'computational noise'. An extra Halley step would double the runtime for a tiny gain and so is not chosen for this implementation, but remains a option, as detailed above.

For the Lambert W_-1 branch, the Fukushima algorithm is used.

Lookup tables

For speed during the bisection, Fukushima's algorithm computes lookup tables of powers of e and z for integral Lambert W. There are 64 elements in these tables. The FORTRAN version (and the C++ translation by Veberic) computed these (once) as static data. This is slower, may cause trouble with multithreading, and is slightly inaccurate because of rounding errors from repeated(64) multiplications.

In this implementation the array values have been computed using Boost.Multiprecision 50 decimal digit and output as C++ arrays 37 decimal digit long double literals using max_digits10 precision

std::cout.precision(std::numeric_limits<cpp_bin_float_quad>::max_digits10);

The arrays are as const and constexpr and static as possible (for the compiler version), using BOOST_STATIC_CONSTEXPR macro. (See lambert_w_lookup_table_generator.cpp The precision was chosen to ensure that if used as long double arrays, then the values output to lambert_w_lookup_table.ipp will be the nearest representable value for the type chose by a typedef in lambert_w.hpp.

typedef double lookup_t; // Type for lookup table (`double` or `float`, or even `long double`?)

This is to allow for future use at higher precision, up to platforms that use 128-bit (hardware or software) for their long double type.

The accuracy of the tables was confirmed using Wolfram Alpha and agrees at the 37th decimal place, so ensuring that the value is exactly read into even 128-bit long double to the nearest representation.

Higher precision

For types more precise than double, Fukushima reported that it was best to use the double estimate as a starting point, followed by refinement using Halley iterations or other methods; our experience confirms this.

Using Boost.Multiprecision it is simple to compute very high precision values of Lambert W at least to thousands of decimal digits over most of the range of z arguments.

For this reason, the lookup tables and bisection are only carried out at low precision, usually double, chosen by the typedef double lookup_t. Unlike the FORTRAN version, the lookup tables of Lambert_W of integral values are precomputed as C++ static arrays of floating-point literals. The default is a typedef setting the type to double. To allow users to vary the precision from float to long double these are computed to 128-bit precision to ensure that even platforms with long double do not lose precision.

The FORTRAN version and translation only permits the z argument to be the largest items in these lookup arrays, wm0s[64] = 3.99049, producing an error message and returning NaN. So 64 is the largest possible value ever returned from the lambert_w0 function. This is far from the std::numeric_limits<>::max() for even floats. Therefore this implementation uses an approximation or 'guess' and Halley's method to refine the result. Logarithmic approximation is discussed at length by R.M.Corless et al. (page 349). Here we use the first two terms of equation 4.19:

T lz = log(z);
T llz = log(lz);
guess = lz - llz + (llz / lz);

This gives a useful precision suitable for Halley refinement.

Similarly, for Lambert W_-1 branch, tiny values very near zero, W > 64 cannot be computed using the lookup table. For this region, an approximation followed by a few (usually 3) Halley refinements. See wm1_near_zero.

For the less well-behaved regions for Lambert W₀ z arguments near zero, and near the branch singularity at -1/e, some series functions are used.

Small values of argument z near zero

When argument z is small and near zero, there is an efficient and accurate series evaluation method available (implemented in lambert_w0_small_z). There is no equivalent for the W_-1 branch as this only covers argument z < -1/e. The cutoff used abs(z) < 0.05 is as found by trial and error by Fukushima.

Coefficients of the inverted series expansion of the Lambert W function around z = 0 are computed following Fukushima using 17 terms of a Taylor series computed using Wolfram Mathematica with

InverseSeries[Series[z Exp[z],{z,0,17}]]

See Tosio Fukushima, Journal of Computational and Applied Mathematics 244 (2013), page 86.

To provide higher precision constants (34 decimal digits) for types larger than long double,

InverseSeries[Series[z Exp[z],{z,0,34}]]

were also computed, but for current hardware it was found that evaluating a double precision and then refining with Halley's method was quicker and more accurate.

Decimal values of specifications for built-in floating-point types below are 21 digits precision == std::numeric_limits<T>::max_digits10 for long double.

Specializations for lambert_w0_small_z are provided for float, double, long double, float128 and for Boost.Multiprecision types.

The tag_type selection is based on the value std::numeric_limits<T>::max_digits10 (and not on the floating-point type T). This distinguishes between long double types that commonly vary between 64 and 80-bits, and also compilers that have a float type using 64 bits and/or long double using 128-bits.

As noted in the implementation section above, it is only possible to ensure the nearest representable value by casting from a higher precision type, computed at very, very much greater cost.

For multiprecision types, first several terms of the series are tabulated and evaluated as a polynomial: (this will save us a bunch of expensive calls to pow). Then our series functor is initialized "as if" it had already reached term 18, enough evaluation of built-in 64-bit double and float (and 80-bit long double) types. Finally the functor is called repeatedly to compute as many additional series terms as necessary to achive the desired precision, set from get_epsilon (or terminated by evaluation_error on reaching the set iteration limit max_series_iterations).

A little more than one decimal digit of precision is gained by each additional series term. This allows computation of Lambert W near zero to at least 1000 decimal digit precision, given sufficient compute time.

Argument z near the singularity at -1/e between branches W₀ and W_-1

Variants of Function lambert_w_singularity_series are used to handle z arguments which are near to the singularity at z = -exp(-1) = -3.6787944 where the branches W₀ and W_-1 join.

T. Fukushima / Journal of Computational and Applied Mathematics 244 (2013) 77-89 describes using Wolfram Mathematica

InverseSeries\[Series\[sqrt\[2(p Exp\[1 + p\] + 1)\], {p,-1, 20}\]\]

to provide his Table 3, page 85.

This implementation used Wolfram Mathematica to obtain 40 series terms at 50 decimal digit precision

N\[InverseSeries\[Series\[Sqrt\[2(p Exp\[1 + p\] + 1)\], { p,-1,40 }\]\], 50\]

-1+p-p^2/3+(11 p^3)/72-(43 p^4)/540+(769 p^5)/17280-(221 p^6)/8505+(680863 p^7)/43545600 ...

These constants are computed at compile time for the full precision for any RealType T using the original rationals from Fukushima Table 3.

Longer decimal digits strings are rationals pre-evaluated using Wolfram Mathematica. Some integer constants overflow, so largest size available is used, suffixed by uLL.

Above the 14th term, the rationals exceed the range of unsigned long long and are replaced by pre-computed decimal values at least 21 digits precision == max_digits10 for long double.

A macro BOOST_MATH_TEST_VALUE (defined in test_value.hpp) taking a decimal floating-point literal was used to allow testing with both built-in floating-point types like double which have contructors taking literal decimal values like 3.14, and also multiprecision and other User-defined Types that only provide full-precision construction from decimal digit strings like "3.14". (Construction of multiprecision types from built-in floating-point types only provides the precision of the built-in type, like double, only 17 decimal digits).

	Tip
	Be exceeding careful not to silently lose precision by constructing multiprecision types from literal decimal types, usually double. Use decimal digit strings like "3.1459" instead. See examples.

Fukushima's implementation used 20 series terms; it was confirmed that using more terms does not usefully increase accuracy.

Lambert W_-1 arguments values very near zero.

The lookup tables of Fukushima have only 64 elements, so that the z argument nearest zero is -1.0264389699511303e-26, that corresponds to a maximum Lambert W_-1 value of 64.0. Fukushima's implementation did not cater for z argument values that are smaller (nearer to zero), but this implementation adds code to accept smaller (but not denormalised) values of z. A crude approximation for these very small values is to take the exponent and multiply by ln[10] ~= 2.3. We also tried the approximation first proposed by Corless et al. using ln(-z), (equation 4.19 page 349) and then tried improving by a 2nd term -ln(ln(-z)), and finally the ratio term -ln(ln(-z))/ln(-z).

For a z very close to z = -1.0264389699511303e-26 when W = 64, when effect of ln(ln(-z) term, and ratio L1/L2 is greatest, the possible 'guesses' are

z = -1.e-26, w = -64.02, guess = -64.0277, ln(-z) = -59.8672, ln(-ln(-z) = 4.0921, llz/lz = -0.0684

whereas at the minimum (unnormalized) z

z = -2.2250e-308, w = -714.9, guess = -714.9687, ln(-z) = -708.3964, ln(-ln(-z) = 6.5630, llz/lz = -0.0092

Although the addition of the 3rd ratio term did not reduce the number of Halley iterations needed, it might allow return of a better low precision estimate without any Halley iterations. For the worst case near w = 64, the error in the 'guess' is 0.008, ratio 0.0001 or 1 in 10,000 digits 10 ~= 4. Two log evalutations are still needed, but is probably over an order of magnitude faster.

Halley's method was then used to refine the estimate of Lambert W_-1 from this guess. Experiments showed that although all approximations reached with Unit in the last place (ULP) of the closest representable value, the computational cost of the log functions was easily paid by far fewer iterations (typically from 8 down to 4 iterations for double or float).

Halley refinement

After obtaining a double approximation, for double, long double and quad 128-bit precision, a single iteration should suffice because Halley iteration should triple the precision with each step (as long as the function is well behaved - and it is), and since we have at least half of the bits correct already, one Halley step is ample to get to 128-bit precision.

Lambert W Derivatives

The derivatives are computed using the formulae in Wikipedia.

Testing

Initial testing of the algorithm was done using a small number of spot tests.

After it was established that the underlying algorithm (including unlimited Halley refinements with a tight terminating criterion) was correct, some tables of Lambert W values were computed using a 100 decimal digit precision Boost.Multiprecision cpp_dec_float_100 type and saved as a C++ program that will initialise arrays of values of z arguments and lambert_W0 (lambert_w_mp_high_values.ipp and lambert_w_mp_low_values.ipp ).

(A few of these pairs were checked against values computed by Wolfram Alpha to try to guard against mistakes; all those tested agreed to the penultimate decimal place, so they can be considered reliable to at least 98 decimal digits precision).

A macro BOOST_MATH_TEST_VALUE was used to allow tests with any real type, both fundamental (built-in) types and Boost.Multiprecision. (This is necessary because fundamental (built-in) types have a constructor from floating-point literals like 3.1459F, 3.1459 or 3.1459L whereas Boost.Multiprecision types may lose precision unless constructed from decimal digits strings like "3.1459").

The 100-decimal digits precision pairs were then used to assess the precision of less-precise types, including Boost.Multiprecision cpp_bin_float_quad and cpp_bin_float_50. static_casting from the high precision types should give the closest representable value of the less-precise type; this is then be used to assess the precision of the Lambert W algorithm.

Tests using confirm that over nearly all the range of z arguments, nearly all estimates are the nearest representable value, a minority are within 1 Unit in the last place (ULP) and only a very few 2 ULP.

For the range of z arguments over the range -0.35 to 0.5, a different algorithm is used, but the same technique of evaluating reference values using a Boost.Multiprecision cpp_dec_float_100 was used. For extremely small z arguments, near zero, and those extremely near the singularity at the branch point, precision can be much lower, as might be expected.

See source at: lambert_w_simple_examples.cpp test_lambert_w.cpp contains routine tests using Boost.Test. lambert_w_errors_graph.cpp generating error graphs.

Testing with quadrature

A further method of testing over a wide range of argument z values was devised by Nick Thompson (cunningly also to test the recently written quadrature routines including Boost.Multiprecision !). These are definite integral formulas involving the W function that are exactly known constants, for example, LambertW0(1/(z²) == √(2π), see Definite Integrals. Some care was needed to avoid overflow and underflow as the integral function must evaluate to a finite result over the entire range.

Other implementations

The Lambert W has also been discussed in a Boost thread.

This also gives link to a prototype version by which also gives complex results (x < -exp(-1), about -0.367879). Balazs Cziraki 2016 Physicist, PhD student at Eotvos Lorand University, ELTE TTK Institute of Physics, Budapest. has also produced a prototype C++ library that can compute the Lambert W function for floating point and complex number types. This is not implemented here but might be completed in the future.

Acknowledgements

References

Initial guesses based on: