Lambert W Function

Source ยป Wikipedia

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\ldots\)
\(W(-1)\)\(-0.31813\ldots+1.33723\ldots i\)
\(W(i)\)\(0.374699\ldots+0.576413\ldots i\)
\(W(-i)\)\(0.374699\ldots-0.576413\ldots i\)
\(W(e)\)\(1\)
\(W(-e)\)\(0.39497\ldots+1.78818\ldots 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}=w_j-\frac{w_je^{w_j}-x}{e^{w_j}+w_je^{w_j}}$$

Alternatively, when simplified: $$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\le10^9\).

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.

Example

The run time was approximately 1 second.

0 <= x <= 1e9
W(x), x = 1

Convergence:
        1
        0.6839397205857211607977618850807304337229055655158839172539184008487307478724499016785736371729598219
        0.5774544771544497158879277307209708744246953720136943334618725831739970553254359545244940341535829055
        0.5672297377301170424199788873305659381989740529912345334258918256943726362413406459666277736815827773
        0.5671432965302959290104061393322652388244468110017561469052367139960666960291476041675058924102131759
        0.5671432904097839036821986273638501915396349920558000709178381348786164164929847042875694697488581049
        0.5671432904097838729999686622103563208086217342023287565768851248326011789240991181028762582526321268
        0.5671432904097838729999686622103555497538157871865125081351310792235327403215041862450323133272582057
        0.5671432904097838729999686622103555497538157871865125081351310792230457930866845666932194469617522946
        0.5671432904097838729999686622103555497538157871865125081351310792230457930866845666932194469617522946

W(1) = 0.5671432904097838729999686622103555497538157871865125081351310792230457930866845666932194469617522946

(rounded up to 100 digits, precise after 9 iterations)

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!

Example

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

Initial Guess: 23.51058492712918919114582243355488790981166285829588560557376358648060096831120133580533938585632116
Convergence:
        23.51058492712918919114582243355488790981166285829588560557376358648060096831120133580533938585632116
        23.51347595210448454142871835682176384958738553442394865049432759042578858344306346044132670280112569
        23.51347595429244717613280522620836780914305001935314712126885835149239110249930323970435009335473405
        23.51347595429244717613280522715678680545105184461531350781886276647715841787449182724224065575351117
        23.51347595429244717613280522715678680545105184461531350781886276647715841787449182731948759374067864

W(382903813902) = 23.51347595429244717613280522715678680545105184461531350781886276647715841787449182731948759374067864
(rounded up to 100 digits, precise after 5 iterations)

Lambert-W-function-v3.cpp

Example

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

Initial Guess: 44.35339134883215996973968716884070132445691959510133480762678367005160538324305405659125047252476026
Convergence:
        44.35339134883215996973968716884070132445691959510133480762678367005160538324305405659125047252476026
        44.3549624063445606412268227676937814316982644538669196388950582583829684964816572864522199228692494
        44.35496240668242513973657868748029822161096853261157563622792382991592856784857054207453437390011557
        44.35496240668242513973657868748365863970197225075861669793254714385770220419310494747348654405936291
        44.35496240668242513973657868748365863970197225075861669793254714385770220419310494747348654736568818

W(812938391290923992913) = 44.35496240668242513973657868748365863970197225075861669793254714385770220419310494747348654736568818
(rounded up to 100 digits, precise after 5 iterations)

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)