Lambert W Function (hp-42s)

[Albert Chan]

A better convergence test might be (x-prev_x)*eps + x - x = 0

With Newton's method, eps can be set to "half" the calculator precision.
Example, 12 digits calculator, eps = 1e-6 is sufficient.

Quote:For example let's consider a infinite tower of the imaginary number i, i^i^i^i......

Instead of using LambertW, we may also do the direct route.

x = k^k^k ... = k^x
log(x) = x log(k)
Let f = log(x) - x log(k), f' = 1/x - log(k)

Newton's iteration x = x - f(x)/f'(x):

x = (1 - log(x)) / (1/x - log(k)
y = 1/x = (y - log(k)) / (log(y) + 1)

>>> from cmath import *
>>> a = -log(1j)                                     # k = 1j
>>> g = lambda y,a: (y+a)/(log(y)+1)     # g(1,a) = 1+a
>>> g(1+a,a) →(1.3127784610672455-1.1245675977983851j)
>>> g(_,a)     →(1.3607126098119278-1.119058042805926j)
>>> g(_,a)     →(1.3606248500898628-1.1194391815927829j)
>>> g(_,a)     →(1.3606248702911174-1.1194391662423497j)
>>> g(_,a)     →(1.3606248702911177-1.11943916624235j), pass convergence test, eps=2**-26
>>> y = _
>>> x = 1/y   →(0.43828293672703211+0.36059247187138549j)
>>> 1j ** x    →(0.43828293672703211+0.36059247187138554j)

x = 1/y = W(a)/a → (a/y) exp(a/y) = a → exp(a/y) = y → W(a) = a/y = log(y)

>>> a/y         → (0.56641733028546448-0.68845322710770207j)
>>> log(y)     → (0.56641733028546437-0.68845322710770218j)

Comment: Solving for x = k^x is not quite the same as x = k^k^k ...
Example:

3 = (³√3)^3, but (³√3) ^ (³√3) ^ (³√3) ... ≈ 2.47805268029 ≠ 3