Lambert W function (for HP Prime)

[lyuka]

(11-05-2020 03:18 PM)Albert Chan Wrote:  Here is proposed patch, for LW0H, Rev 1.36 (Nov 5, 2020):

Code:

    ...
    y := y - t / v;  // Householder's method of order 4 
    t := abs(y - p); // correction radius
  UNTIL 0 == t OR (r < 1 AND r <= t); // convergence check
  return p;
END

Thank you very much for your analysis and your patch proposal.

The code updated.

Code:

//LW0H - Principal branch of the Lambert W function
//Rev.1.37 (Nov. 6, 2020) (c) Takayuki HOSODA (aka Lyuka)
//Acknowledgments: Thanks to Albert Chan for his informative suggestions.
EXPORT LW0H(x)
BEGIN
  LOCAL y, p, r, s, t, u, v;
  HComplex := 1; // Enable complex result from a real input.
  r := 1 / e;
  IF (x == -r) THEN return -1; END;
  s := x + r;
  t := e - sqrt(2) - 1;
  y := ln(r + sqrt((r + r) * s) + t * s); // approximation near x=-1/e
  t := 9;
  REPEAT
    p := y;
    r := t;
    s := exp(y);
    u := y * s;
    v := u - x;
    t := s + x;
    v := 6 * ((y + 1) * u * x + t * (t + u)) + v * (3 * (u + x) + y * v);
    IF (0 == v) THEN break; END;
    t := 3 * (u - x) * ((u + x) * (y + 2) + s + s);
    y := y - t / v; // Householder's method of order 4
    t := abs(y - p); // correction radius;
  UNTIL 0 == t OR (r < 1 AND r <= t); // convergence check
  return p;
END;

LW0H Rev.1.37 convergence example at x = -1.78