Lambert W Function (hp-42s)

[Albert Chan]

Hi, Juan14

I was experimenting another way to get W(a), via infinite tetration route.
Solving for y = exp(W(a)) maybe better than directly going for W(a):

Code:

from cmath import *
g = lambda y,a: (y+a)/(log(y)+1)
w = lambda x,a: x+(a/exp(x)-x)/(x+1)

def fixedpoint(f, x, a, eps=2**-26, limit=100):
    for t in xrange(1,limit):   # f() calls
        prev, x = x, f(x,a)
        if (x-prev)*eps + x == x: return x, t

def test(a):
    y, t = fixedpoint(g, 1+a, a)
    print 'y->w =', (a/y, t)
    print 'w    =', fixedpoint(w, log(1+a), a)

>>> test(1)
y->w = ((0.56714329040978384+0j), 4)
w      = ((0.56714329040978384+0j), 5)
>>> test(1e10)
y->w = ((20.028685413304952+0j), 5)
w      = ((20.028685413304952+0j), 8)
>>> test(1e100)
y->w = ((224.84310644511851+0j), 4)
w      = ((224.84310644511851+0j), 10)

Note: both versions use the same guess, y = e^x = 1+a
Note: y->w uses a/y instead of log(y) to recover W(a), to reduce errors when y ≈ 1