for HP Prime
An implementation of
Lambert W function and exponent of Lambert W function
for HP Prime

Oct. 25 2020
Takayuki HOSODA

Lambert W function exponent of Lambert W function

From the relation below, the Lambert W can be calculated by using eW — the exponent of Lambert W function.

W(x) = \frac{x}{\mathrm{e}^{W(x)}} = \ln\big(\mathrm{e}^{W(x)}}\big\,)
eW0 — Exponent of the principal branch of the Lambert W function

eW0(x) — the exponent of the principal branch of the Lambert W function W0(x) — can be calculated by using Newton-Raphson method,

y_\mathrm{n+1} \leftarrow y_\mathrm{n} - \frac{y_\mathrm{n}\ln(y_\mathrm{n}) - x}{\ln(y_\mathrm{n}) + 1}
or Halley's method shown below.
y_\mathrm{n+1} \leftarrow y_\mathrm{n} + \frac {2y_\mathrm{n}(1 + \ln y_\mathrm{n}) \,(x - y_\mathrm{n} \ln y_\mathrm{n})} {2y_\mathrm{n}(1 + \ln y_\mathrm{n})^2 + x - y_\mathrm{n} \ln y_\mathrm{n}}

When finding eW using an iterative method such as Newton's Rapson method, giving a bad initial value, especially when it is very close to the singularity x = -1 / e, can be terrible in terms of convergence. Therefore, it is desirable that the initial value asymptotically approaches eW at the singularity. For that reason, I chose the following formula as an approximate expression that gives the initial value. This equation replaces the coefficient 1 / 3 of the third term of the Puiseux series expansion of e W0 at x = -1 / e with 0.3 and ignores the fourth and subsequent terms.

Initial approximation :

y_0 = \frac{1}{\mathrm{e}} + \sqrt{\frac{2}{\mathrm{e}} \, \bigg(\frac{1}{\mathrm{e}} + x\bigg)} + 0.3 \,\bigg(\frac{1}{\mathrm{e}} + x\bigg)
Series expansion of eW at z = -1 / e (Puiseux series) :
\frac{1}{\mathrm{e}} + \sqrt{\frac{2}{\mathrm{e}}\left(z + \frac{1}{\mathrm{e}}\right)}} +\frac{1}{3} \left(z + \frac{1}{\mathrm{e}}\right) - \frac{1}{18} \sqrt{\frac{\mathrm{e}}{2}} \left(z + \frac{1}{\mathrm{e}}\right)^{\big{\frac{3}{2}}} +O(z^2)
Instead of the coefficient of 0.3 in the formula above, it can be (e - √2 - 1) which make it zero at x = 0,
y_0 = \frac{1}{\mathrm{e}} + \sqrt{\frac{2}{\mathrm{e}} \, \bigg(\frac{1}{\mathrm{e}} + x\bigg)} + (\mathrm{e} - \sqrt{2} - 1) \,\bigg(\frac{1}{\mathrm{e}} + x\bigg)
or "1 / 3" from Puiseux series, but 0.3 is quite good enough for this purpose.

It's nice to be able to handle complex numbers easily with HP Prime, and thanks to sqrt √ in this equation, HP Prime can automatically get the complex y0 from x less than -1 / e.
Note that, complex output from real input must be allowed by selecting the check box of 'Complex' in 'Home Setting'.

eW (for HP Prime)
eW(x) returns the requested exponent of the principal branch of the Lambert W function.

Recurrence formula (Halleys method) :

    y_\mathrm{n+1} \leftarrow y_\mathrm{n} + \frac
{2y_\mathrm{n}(1 + \ln y_\mathrm{n}) \,(x - y_\mathrm{n} \ln y_\mathrm{n})}
{2y_\mathrm{n}(1 + \ln y_\mathrm{n})^2 + x - y_\mathrm{n} \ln y_\mathrm{n}}

Initial approximation :

    y_0 = \frac{1}{\mathrm{e}} + \sqrt{\frac{2}{\mathrm{e}} \, \bigg(\frac{1}{\mathrm{e}} + x\bigg)}  + 0.3 \,\bigg(\frac{1}{\mathrm{e}} + x\bigg)

    

#CAS
//eW : exponent of the Lambert W0 function
//Rev.1.52 (Dec. 4, 2020) (c) Takayuki HOSODA (aka Lyuka)
// http://www.finetune.co.jp/~lyuka/technote/lambertw/
//Acknowledgments: Thanks to Albert Chan for his informative suggestions. 
eW(x):=
BEGIN
  LOCAL y, p, s, t, u, v;
  LOCAL q, r, m;
  HComplex := 1; // Enable complex result from a real input.
  r := 1.0 / e;
  q := e - sqrt(2) - 1;
  s := x + r;
  IF s == 0 THEN return r; END;
  m := MAXREAL / 1024.0;
  IF abs(RE(x)) <= m AND abs(IM(x)) <= m THEN
    m := 1.0;
    y := r + sqrt(2.0 * r * s) + q * s; // approximation near x=-1/e
  ELSE
    m := 1.0 / 1024.0;
    x := x * m;
    s := x * m;
    y := r * m + sqrt(r * (s * 64.0)) + q * s;
    y := y * 1024.0;
  END;
  r := MAXREAL;
  REPEAT
    q := r;
    p := y;
    t := ln(y);
    v := 1 + t;
    s := 0.5 * (x * inv(y) - t * m) * inv(v) + v * m; //s <> 0, y <> 0, v <> 0
    t := (x - y * (t * m));
    y := y + t * inv(s); // Halley's method
    r := cabs(y - p); // correction radius
  UNTIL 0 == r OR q <= r; // convergence check
  return p;
END;
// This is a workaround for the dynamic range problem of complex division and abs() in CAS.
// Conforms to HP Prime software version 2.1.14433 (2020 01 21).
// cabs(z) -- returns the magnitude of the complex number z.
cabs(z):=
BEGIN
  LOCAL x, y, t, s;
  x := RE(z);
  y := IM(z);
  IF x < 0 THEN x := -x; END;
  IF y < 0 THEN y := -y; END;
  IF x == 0 THEN return y; END;
  IF y == 0 THEN return x; END;
  IF (y > x) THEN
    t := x; x := y; y := t;
  END;
  t := x - y;
  IF t > y THEN
    t := x / y;
    t := t + sqrt(1 + t * t);
  ELSE
    t := t / y;
    s := t * (t + 2);
    t := t + s / (sqrt(2) + sqrt(2 + s));
    t := t + sqrt(2) + 1;
  END;
  return(x + y / t);
END;
#END
Download 'eW.hpprgm' for HP Prime.
Calculation results by WolframAlpha
W0 — The principal branch of the Lambert W function

It is necessary to consider the singularity, but since eW does not become zero, it was relatively easy to find eW by Newton-Raphson method or Harrey's method. However, when trying to find W0 - Lambert's W of the primary branch by such an iterative method, since W0 becomes zero when x is zero, so it is necessary to consider zero as well as the singularity. In other words, it is strongly desired that the approximation formula of the initial value of the iterative method has an asymptotic to zero error at the singular point and zero. Also, since W0 is, from the singularity, a monotonically increasing function, its approximation must also be a monotonically increasing function. To meet these requirements, the initial value is calculated using the following approximation function.

y_0 = \frac{1}{\mathrm{e}} + \sqrt{\frac{2}{\mathrm{e}} \, \bigg(\frac{1}{\mathrm{e}} + x\bigg)} + (\mathrm{e} - \sqrt{2} - 1) \,\bigg(\frac{1}{\mathrm{e}} + x\bigg)

Since W0 has a large derivative around the singularity, Newton's method can cause catastrophic behavior near the singularity if the initial approximation is poor or due to computational errors. Therefore, we will use the Halley method or Householder's method of order 4, which is more stable than the Newton method, to obtain W0.

In addition to the Halley's method of cubic convergence, the Householder's method of order 4 can also be used to calculate W0. However, special care is required because the calculation is complicated and it is easy for digits to overflow during the calculation. Householder's method was only a few percent faster than Halley's method for calculating up to 16 decimal places.

 Newton-Raphson method :
y_\mathrm{n+1} \leftarrow y_\mathrm{n} + \frac {x - y_\mathrm{n}\, \mathrm{e}^{\, \displaystyle y_\mathrm{n}}} { \mathrm{e}^{\,\displaystyle y_\mathrm{n}} + y_\mathrm{n}\,\mathrm{e}^{\, \displaystyle y_\mathrm{n}}
 Halley's method :
y_\mathrm{n+1} \leftarrow y_\mathrm{n} - \frac {2 (y_\mathrm{n} + 1) \, (\mathrm{e}^{\, \displaystyle y_\mathrm{n}} - x)} { ({ y_\mathrm{n}}^2 + 2 y_\mathrm{n} + 2) \, \mathrm{e}^{\,\displaystyle y_\mathrm{n}} + (y_\mathrm{n} + 2) x}
 Householder's method of order 4 :
y_{\mathrm n+1} \leftarrow y_{\mathrm n} - \frac{3y_{\mathrm n}^3(\mathrm{e}^{\, \displaystyle y_\mathrm{n}})^2+6y_{\mathrm n}^2(\mathrm{e}^{\, \displaystyle y_\mathrm{n}})^2-3y_{\mathrm n}{x}^2+6y_{\mathrm n}(\mathrm{e}^{\, \displaystyle y_\mathrm{n}})^2-6{x}^2-6{x}\mathrm{e}^{\, \displaystyle y_\mathrm{n}}} {{6y_{\mathrm n}}^2{x}\mathrm{e}^{\, \displaystyle y_\mathrm{n}}+{y_{\mathrm n}}^2\mathrm{e}^{\, \displaystyle y_\mathrm{n}}+18y_{\mathrm n}{x}\mathrm{e}^{\, \displaystyle y_\mathrm{n}}+6y_{\mathrm n}(\mathrm{e}^{\, \displaystyle y_\mathrm{n}})^2+5y_{\mathrm n}\mathrm{e}^{\, \displaystyle y_\mathrm{n}}+12{x}\mathrm{e}^{\, \displaystyle y_\mathrm{n}}+6(\mathrm{e}^{\, \displaystyle y_\mathrm{n}})^2+6\mathrm{e}^{\, \displaystyle y_\mathrm{n}}}

The convergence test by Urabe's theorem is used to determine the convergence. It is a method to judge that the calculation limit has been reached and terminate the iteration when the correction value does not change or increases.

W0 calculation using Halley's method
LW0C (for HP Prime CAS)
LW0C(x) returns the requested value of the principal branch of the Lambert W function.

Recurrence formula (Halley's method) :

    y_\mathrm{n+1} \leftarrow y_\mathrm{n} - \frac
{2 (y_\mathrm{n} + 1) \, (\mathrm{e}^{\, \displaystyle y_\mathrm{n}} - x)}
{ ({ y_\mathrm{n}}^2 + 2 y_\mathrm{n} + 2) \, 
\mathrm{e}^{\,\displaystyle y_\mathrm{n}} + (y_\mathrm{n} + 2) x}

Initial approximation :

    y_0 = \ln\left(\frac{1}{\mathrm{e}} + \sqrt{\frac{2}{\mathrm{e}} \, \left(\frac{1}{\mathrm{e}} + x\right)}  + (\mathrm{e} -  \sqrt{2} - 1) \,\left(\frac{1}{\mathrm{e}} + x\right)\right)

    


#CAS
//LW0C - Principal branch of the Lambert W function
//Rev.1.52 (Dec. 4, 2020) (c) Takayuki HOSODA (aka Lyuka)
//Acknowledgments: Thanks to Albert Chan for his informative suggestions.
LW0C(x):=
BEGIN
  LOCAL y, p, s, t;
  LOCAL q, r, m;
  HComplex := 1; // Enable complex result from a real input.
  r := 1.0 / e;
  q := e - sqrt(2.0) - 1.0;
  m := MAXREAL / 1024.0;
  IF abs(RE(x)) <= m AND abs(IM(x)) <= m THEN
    m := 1.0;
    s := x + r;
    IF s == 0 THEN return -1; END;
    y := ln(r + sqrt(2.0 * r * s) + q * s); // approximation near x=-1/e
  ELSE
    m := 1.0 / 1024.0;
    x := x * m;
    s := x * m;
    y := ln(r * m + sqrt(r * (s * 64.0)) + q * s);
    y := y + 6.931471805599453094; // + ln(2) * 10
  END;
  r := MAXREAL; //1.79769313486e308
  REPEAT
    q := r;
    p := y;
    t := exp(y) * m;
    s := y * t - x;
    t := (0.5 * y) * inv(y + 1.0) * s + t  + x; // t <> 0
    y := y - s * inv(t); // Halley's method
    r := abs(y - p); // correction radius;
  UNTIL 0 == r OR q <= r; // convergence check
  return p;
END;
#END
Download 'LW0C.hpprgm' for HP Prime.

Calculation example of LW0C

Calculation results by WolframAlpha W0 calculation using Householder's method of order 4
W0 by Using Householder's method of order 4
LW0H (for HP Prime)
LW0H(x) returns the requested value of the principal branch of the Lambert W function.

Recurrence formula (Householder's method of order 4) :

   y_{\mathrm n+1} \leftarrow y_{\mathrm n} - 
\frac{3y_{\mathrm n}^3(\mathrm{e}^{\, \displaystyle y_\mathrm{n}})^2+6y_{\mathrm n}^2(\mathrm{e}^{\, \displaystyle y_\mathrm{n}})^2-3y_{\mathrm n}{x}^2+6y_{\mathrm n}(\mathrm{e}^{\, \displaystyle y_\mathrm{n}})^2-6{x}^2-6{x}\mathrm{e}^{\, \displaystyle y_\mathrm{n}}}
{{6y_{\mathrm n}}^2{x}\mathrm{e}^{\, \displaystyle y_\mathrm{n}}+{y_{\mathrm n}}^2\mathrm{e}^{\, \displaystyle y_\mathrm{n}}+18y_{\mathrm n}{x}\mathrm{e}^{\, \displaystyle y_\mathrm{n}}+6y_{\mathrm n}(\mathrm{e}^{\, \displaystyle y_\mathrm{n}})^2+5y_{\mathrm n}\mathrm{e}^{\, \displaystyle y_\mathrm{n}}+12{x}\mathrm{e}^{\, \displaystyle y_\mathrm{n}}+6(\mathrm{e}^{\, \displaystyle y_\mathrm{n}})^2+6\mathrm{e}^{\, \displaystyle y_\mathrm{n}}}


Initial approximation :

         y_0 = \ln\left(\frac{1}{\mathrm{e}} + \sqrt{\frac{2}{\mathrm{e}} \, \left(\frac{1}{\mathrm{e}} + x\right)}  + (\mathrm{e} -  \sqrt{2} - 1) \,\left(\frac{1}{\mathrm{e}} + x\right)\right)


#CAS
//LW0H - Principal branch of the Lambert W function
//Rev.1.46 (Nov. 28, 2020) (c) Takayuki HOSODA (aka Lyuka)
//Acknowledgments: Thanks to Albert Chan for his informative suggestions.
LW0H(x):=
BEGIN
  LOCAL y, p, s, t, u, v;
  LOCAL q, r;
  HComplex := 1; // Enable complex result from a real input.
  r := sqrt(MAXREAL) / 1e4; // Limit x range
  IF abs(RE(x)) > r OR abs(IM(x)) > r THEN return "LW0H: x out of range"; END;
  r := 1.0 / e;
  s := x + r;
  IF s == 0 THEN return -1; END;
  q := e - sqrt(2.0) - 1.0;
  y := 1.0 * ln(r + sqrt(2.0 * r * s) + q * s); // approximation near x=-1/e
  r := MAXREAL;
  REPEAT
    q := r;
    p := y;
    s := exp(y);
    t := s + x;
    u := y * s;
    v := u - x;
    v := (6.0 * ((y + 1.0) * u * x + t * (t + u)) + v * (3.0 * (u + x) + y * v));
    IF (0 == v) THEN break; END;
    t := 3.0 * (u - x) * ((u + x) * (y + 2.0) + 2.0 * s);
    y := y - t * inv(v); // Householder's method of order 4
    r := abs(y - p); // correction radius
  UNTIL 0 == r OR q <= r; // convergence check by Urabe's theorem
  return p;
END;
#END
Download 'LW0H.hpprgm' for HP Prime.

Calculation speed comparison, eW (Rev.1.46) V.S. LW0C (Rev.1.46) V.S. LW0H (Rev.1.46)

Download 'scLW.hpprgm' used to make the comparison above -- invoked as scLW(3) or scLW(50).
From some screen captures of HP Prime for reference



Reference
  1. "Lambert W-function", Wolfram MathWorld
  2. "Why W?", Hayes, B. (2005). American Scientist. 93 (2): 104‐108.
  3. "Lambert W function", Wikipedia
  4. "The Lambert W function poster", Ontario Research Centre for Computer Algebra.
  5. "Puiseux series", Wikipedia
  6. "Newton's Method", Wolfram MathWorld
  7. "Householder's Method", Wolfram MathWorld
  8. "Convergence of Numerical Iteration in Solution of Equations.", Urabe Minoru, J. Sci. Hiroshima Univ. Ser. A 19 (1956), no. 3, 479--489. doi:10.32917/hmj/1556071264.
  9. "Improved Complex Division", Michael Baudin, Robert Smith, 2011
  10. GSL - GNU Scientific Library
  11. gnuplot demo - 3D plotting complex functions with color mapping using pm3d
  12. "Mathematics Written in Sand", W. Kahan, University of California @ Berkeley, Nov. 1983
External Links See Also
www.finetune.co.jp [Mail] © 2000 Takayuki HOSODA.