HP49-50G Lambert function

[Albert Chan]

(10-19-2020 09:13 AM)Gil Wrote:  The program gives almost instantaneous answer for non-complex solutions ; as expected,
however, longer times are required for the iterations when dealing with complex numbers solutions.

You can speed up convergence by solving y = e^W(x) instead.

W e^W = y * log(y) = x

f = y * log(y) - x, f' = log(y) + 1, then apply Newton's method.

f' is insensitive to guess y, and will converge very fast, even with bad guess.

Code:

from cmath import *
def eW(x, n=5, r=1/e, y=None):
    if not y: y = 0.3*(x+r) + sqrt(2*r*(x+r)) + r
    if y == r: y += 1e-10   # avoid 0 slope
    for i in range(n): print i,y; y = (y+x)/(log(y)+1)
    return y

>>> eW(-1+2j)
0 (0.912470160638+1.60208654199j)
1 (0.985261305638+1.5911210409j)
2 (0.985506513217+1.59184289853j)
3 (0.98550646367+1.59184283455j)
4 (0.98550646367+1.59184283455j)
(0.98550646366951244+1.5918428345475875j)

>>> eW(-1e10+2e10j)
0 (-2999932566.3+6000109109.25j)
1 (-452114322.526+1139389284.48j)
2 (-387942931.588+998560851.502j)
3 (-387765196.033+998129458.041j)
4 (-387765194.562+998129453.635j)
(-387765194.56226903+998129453.63483131j)

https://www.hpmuseum.org/forum/thread-15...#pid136892
http://www.finetune.co.jp/~lyuka/technot...w-42s.html