Lambert W function (for HP Prime)

[toml_12953]

Here's an ANSI/ISO BASIC adaptation of a C program from

https://stackoverflow.com/questions/6021...in-c-sharp

Code:

DECLARE EXTERNAL FUNCTION LambertW

PRINT LambertW(MAXNUM)

END 

EXTERNAL FUNCTION LambertW(x)
! LambertW is not defined in this section
IF (x < EXP(-1)) THEN
   PRINT "The LambertW-function is not defined for";x;"."
   EXIT FUNCTION
END IF

! computes the first branch for real values only

! amount of iterations (empirically found)
LET amountOfIterations = MAX(4, CEIL(LOG10(x) / 3))

! initial guess is based on 0 < ln(a) < 3
LET w = 3 * LOG(x + 1) / 4 

! Halley's method via eqn (5.9) in Corless et al (1996)
FOR i = 0 TO amountOfIterations-1
   LET w = w - (w * EXP(w) - x) / (EXP(w) * (w + 1) - (w + 2) * (w * EXP(w) - x) / (2 * w + 2))
NEXT i

LET LambertW = w
END FUNCTION