Lambert W function (for HP Prime)

[Albert Chan]

(10-27-2020 04:01 PM)Albert Chan Wrote:  For W, Halley's method is safer than Newton's method. (safest is e^W approach)

We can get the safety of e^W approach, without recovering W at the end.
Since w_guess(0) = W(0), we can remove the case for W(a = 0):

x * exp(x) = a,       a ≠ 0
ln(x) + x = ln(a)

f = x + ln(x/a) = 0
f' = 1 + 1/x = (x+1)/x
f'' = -1/x^2
f''/f' = -1/(x*x+x)

This is version 2, using above f setup:

Code:

w_newton2_corr = lambda x, a: (x + log(x/a)) * x/(x+1)
w_newton2 = lambda x, a: x - w_newton2_corr(x,a)
w_halley2 = lambda x, a: x - 1/(1/w_newton2_corr(x,a) + .5/(x*x+x))

def nest(func, x, a, n=6):
    for i in range(n):
        print i,x;
        try: x = func(x,a)
        except ZeroDivisionError: break

Now, redo the bad guess example, and see the improvements:

>>> from mpmath import *
>>> x = -1e99
>>> w = 200 + 3j # bad guess

>>> nest(w_newton2, w, x)
0 (200+3j)
1 (222.54478620735 + 3.12764616976661j)
2 (222.550670661156 + 3.12754041757397j)
3 (222.550670661503 + 3.12754041755172j) <-- converged in 3
4 (222.550670661503 + 3.12754041755172j)
5 (222.550670661503 + 3.12754041755172j)

>>> nest(w_halley2, w, x)
0 (200+3j)
1 (222.551107406185 + 3.12752854221077j)
2 (222.550670661503 + 3.12754041755172j) <-- converged in 2
3 (222.550670661503 + 3.12754041755172j)