Lambert W function (for HP Prime)

[lyuka]

Hi guys,
a few days ago my wife bought me an HP Prime (2AP18AA) as a little early birthday present.
I'm completely new to HP Prime, but I implemented the Lambert W function and the expornent of the Lambert W function that I wrote the other day as a trial to get started.

Please refer to the page below for details.
Lambert W function (for HP Prime)

Code:


//W0 - Principal branch of the Lambert W function
//Rev.1.1 (Oct. 25, 2020) (c) Takayuki HOSODA (aka Lyuka)
EXPORT W0(x)
BEGIN
  LOCAL y, r, t, i;
  r := 1 / e;
  IF (x == -r) THEN
    return -1;
  ELSE
    t := x + r;
    y := ln(r + sqrt((r + r) * t) + 0.3040682660859502 * t); // approximation near x=-1/e
    FOR i FROM 0 TO 16 DO
      r := t;
      t := exp(y);
      IF (-1 == y) THEN break; END;
      t := (x - y * t) / (t * (1.0 + y)); // Newton-Raphson method
      y := y + t;
      IF (abs(t) >= abs(r) AND i > 0) THEN break; END; // convergence check
    END;
  END;
  return y;
END;

Code:


//eW : exponent of the Lambert W0 function
//Rev.1.1 (Oct. 25, 2020) (c) Takayuki HOSODA (aka Lyuka)
EXPORT eW(x)
BEGIN
  LOCAL y, r, t, u, v, i;
  r := 1 / e;
  IF (x == -r) THEN
    return r;
  ELSE
    t := x + r;
    y := r + sqrt((r + r) * t) + 0.3040682660859502 * t; // approximation near x=-1/e
    FOR i FROM 0 to 10 DO
      r := t;
      t := ln(y);
      v := x - y * t;
      t := t + 1;
      u := (t + t) * y;
      t := v + u * t;
      IF (0 == t) THEN break; END;
      t := u * v / t; // Halley's method
      y := y + t;
      IF (abs(t) >= abs(r) AND i > 0) THEN break; END; // convergence check
    END;
  END;
  return y;
END;