Lambert W Function (hp-42s)

[Albert Chan]

(09-28-2020 04:06 PM)lyuka Wrote:  It's nice to be able to handle complex numbers easily with 42S,
but if you try to find e^W with Newton-Raphson method, it will fail very close to -1/e.

I use e^W(x) guess of LN(1+x), it seems OK with close to -1/e

Werner improved on it by doing everything in the stack

Bonus: resulting stack X = Y = eW, Z = T = x. To recover W, we can do either "LN", or "R↓ ÷"

With Free42, tried x = -1/e:

1 [+/-] [e^X] [+/-] XEQ "eW"       → 0.367879441171

It worked, but result "only" have 17 good digits.

Trying it in Python, eW convergence around -1/e is bad, even with good guess.

>>> from cmath import *
>>> g = lambda y,a: (y+a)/(log(y)+1)
>>> x = -1/e
>>> y = log(1+x) # eW guess
>>> y
(-0.45867514538708193+0j)

>>> for i in range(1,101): y = g(y,x); print i, y
...
1 (-0.0183829712865+0.261809735996j)
2 (0.19954449329+0.194333174479j)
3 (0.292276297365+0.112700062519j)
4 (0.332518371694+0.0601902170855j)
5 (0.350846602748+0.0310495664127j)
6 (0.359529649069+0.0157626859032j)
7 (0.363746790135+0.00794072878334j)
8 (0.365823750457+0.00398519956333j)
9 (0.36685426371+0.00199630716024j)
...
100 (0.367879435928+2.54426772287e-11j)