Lambert W function (for HP Prime)

[lyuka]

(10-25-2020 07:02 PM)Albert Chan Wrote:  Does HP Prime have LambertW built-in ?

XCas> solve(x^2 = 2^x, x)

[-exp(-LambertW(ln(2)/2)), exp(-LambertW(-ln(2)/2)), exp(-LambertW(-ln(2)/2,-1))]

XCas> float(ans())

[-0.766664695962, 2.0, 4.0]

---

I think it is safe to use Halley's method for W.
Think of Halley's method as correction to Newton's method slope.
For W, slope correction is at most 50% (when x is *far* from -1/e)

HP Prime doesn't seem to have LambertW function. I found a previous thread on this.
Xcas/Giac vs HP Prime: LambertW

--
I reconfirmed the Halley method.
The precision of the least significant digit is often better than Newton's method,
and I haven't encountered any unstable behavior so far,
so I changed the program to use Halley's method.
The program name has been changed from W0 to LW0.

Code:


//LW0 - Principal branch of the Lambert W function
//Rev.1.22 (Oct. 30, 2020) (c) Takayuki HOSODA (aka Lyuka)
EXPORT LW0(x)
BEGIN
  LOCAL y, r, s, t, i;
  r := 1 / e;
  IF (x == -r) THEN return -1; END;
  t := x + r;
  y := ln(r + sqrt((r + r) * t) + 0.3040682660859502 * t); // approximation near x=-1/e
  FOR i FROM 0 TO 10 DO
    r := t;
    t := exp(y);
    s := y * t;
    t := (1 + 0.5 * y) * (s + x) + t;
    IF (t == 0) THEN break; END;
    t := (y + 1) * (s - x) / t; // Halley's method
    y := y - t;
    IF (t == 0 OR abs(t) >= abs(r) AND i > 0) THEN break; END; // convergence check
  END;
  return y;
END;

Please refer to the page below for details.
Lambert W function and exponent of Lambert W function for HP Prime