Esote esote.net

Lambert W function

The Lambert W function W(x), also called the omega function or the product logarithm, is the inverse of the function f(W)=We^W. It is used as the solution to many problems that cannot be solved using normal methods. While the Lambert W function is actually a set of functions, this program focuses on its principal branch W_0.

Example values of the Lambert W function
W(x) Result
W(0)0
W(1)0.56714
W(-1)-0.31813+1.33723 i
W(i)0.374699+0.576413 i
W(-i)0.374699-0.576413 i
W(e)1
W(-e)0.39497...+1.78818...i

The Lambert W function can be approximated using Newton's method and Halley's method, although methods with faster approximation can be constructed. These methods use successive approximation to converge towards the value of W(x). Newton's method converges quadratically and Halley's method converges cubically.

I do not include imaginary results. I take advantage of Boost to allow for arbitrary precision. You can change the value of const int PRECISION to increase the output precision. Often the last digit is rounded or not precise. This inaccuracy can be counteracted by increasing PRECISION.

Lambert-W-function-v1.cpp

This first version uses Newton's method to converge towards W(x) using the formula:

w_{j+1}=\frac{xe^{-w_j}+w_j^2}{w_j+1}

Breakdown:

  1. Gets the value of x in the form of a string;

  2. Checks if the string contains more than one decimal mark;

  3. Checks if the string is in the format of a decimal fraction (so not the letter "j", for example) by converting it to a double in a try-catch block;

  4. Converts it to the arbFloat data type, the name I chose for Boost's multiprecision cpp_dec_float;

  5. Number is checked to see if it is less than zero (which would return a complex number);

  6. Calculations (explained a few lines down).

In this version, there is a bound on the input of x ≤ 1E9. Inside of the for loop, it is simple:

  1. Convert initial value of wnew to a string, and set the string's precision;

  2. Use the simplified Newton's method to approximate W(x);

  3. Convert the value of wnew after calculations to a string, and set the strings precision;

  4. Compare the pre-calculation string and post-calculation string: if they are the same, it means the precision was exceeded.

Lambert-W-function-v2.cpp

With the new improvements, the value of W(x) within a computable range can be calculated in seconds with around 5 iterations!

Lambert-W-function-v3.cpp

Lambert-W-function-v4.cpp

Example

0 <= x < inf
W(x), x = 58491.4859194

Initial Guess: 8.79912818825742701263554005635953367481999598927155403916723877200983598181966539704818806245543085
Convergence:
	8.79912818825742701263554005635953367481999598927155403916723877200983598181966539704818806245543085
	8.801692453591711067034391690021038727530168930812754421891411409540457859585438570795243845704890877
	8.801692455327469805597939732638277031925622767561811681524957334633929913212419954962801678030124737
	8.801692455327469805597939733176608343067236688470654778832038917101125243949975035911147451410329376
	8.801692455327469805597939733176608343067236688470654778832038917101125243949975035927206894516430969

W(58491.4859194) = 8.801692455327469805597939733176608343067236688470654778832038917101125243949975035927206894516430969
(rounded up to 100 digits, precise after 5 iterations)