Lambert W function (for HP Prime)

[Albert Chan]

(10-27-2020 03:43 AM)lyuka Wrote:  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.

Good call.
For W, Halley's method is safer than Newton's method. (safest is e^W approach)

Newton's method: x ← x - N, where N = f/f'
For W(a), we setup f = x*e^x - a, f' = x*e^x + e^x

XCas> N := (x*e^x - a) / (x*e^x + e^x)
XCas> limit(N, x=+inf)   → 1
XCas> limit(N, x=-inf)    → inf
XCas> limit(N, x=-1)      → inf

Halley's method: x ← x - H,   where H = f / (f' - (f/f')/2*f'') = N / (1 - N/2*(f''/f'))
f'' = (f')' = (f + e^x - a)' = f' + e^x
f''/f' = 1 + e^x/f' = 1 + 1/(1+x)

XCas> H := N/(1 - N/2*(1+1/(1+x)))
XCas> limit(H, x=+inf)   → 2
XCas> limit(H, x=-inf)    → -2
XCas> limit(H, x=-1)      → 0

With bad guess, Newton's correction close to singularity may hugely overshoot.
And, it does not take much to have a bad guess, see the tests below ...

Code:

from mpmath import *
r, c = 1/e, e-sqrt(2)-1

def w_newton_corr(x, a):
    ex = exp(x)
    return (x*ex - a) / (x*ex + ex)

w_guess = lambda x,c=c: log(sqrt(2*r*(x+r)) + c*(x+r) + r)
w_newton = lambda x, a: x - w_newton_corr(x,a)
w_halley = lambda x, a: x - 1/(1/w_newton_corr(x,a) - (.5+.5/(1+x)))

>>> def nest(func, x, a, n=6):
...          for i in range(n): print i,x; x = func(x,a)

>>> x = -1e99
>>> w = 200+3j # bad guess, but only 10% off, W(x) = 222.550670661503 + 3.12754041755172j

>>> nest(w_newton , w, x) # Newton totally failed
0 (200+3j)
1 (6829135788.03541 + 869691954.442808j)
2 (6829135787.03541 + 869691954.442808j)
3 (6829135786.03541 + 869691954.442808j)
4 (6829135785.03541 + 869691954.442808j)
5 (6829135784.03541 + 869691954.442808j)

>>> nest(w_halley, w, x) # 14 iterations to converge
0 (200+3j)
1 (201.990101192677 + 3.00014701205407j)
2 (203.98029891143 + 3.00029117778268j)
3 (205.970591273429 + 3.00043258190625j)
4 (207.960976300103 + 3.00057132498455j)
5 (209.951450963988 + 3.00070764455423j)

Same behavior (but no complex numbers) if we try x = 1e99, guess w = 200.

---

>>> x = -1/e + 1e-5 # close to singularity
>>> w = -1 + 1e-10 # guess too close to singularity, W(x) = -0.992644755197123

>>> nest(w_newton, w, x) # 272,000 iterations to converge
0 -0.9999999999
1 271827.05391625
2 271826.053919929
3 271825.053923608
4 271824.053927287
5 271823.053930966

>>> nest(w_halley, w, x) # 19 iterations to converge
0 -0.9999999999
1 -0.9999999997
2 -0.9999999991
3 -0.9999999973
4 -0.999999991899999
5 -0.999999975699998

With good guess, w_guess(x) = -0.992645539358962, we don't have above problems.
w_newton will converge in 2 iterations, w_halley in 1.