Lambert W Function (hp-42s)

[Albert Chan]

Here is a demo of previous post results, r + err_r = 1/e (108-bits accurate)

Code:

r, err_r = 0.36787944117144233, -1.2428753672788363E-17
g0 = lambda y,a: (y+a)/(log(y)+1)
g1 = lambda y,a: (y+a)/log1p(e*(y - r))
g2 = lambda y,a: (y+a)/log1p(e*(y - r - err_r))

def eW(x, g, n=5):
    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 = g(y,x)
    return y

>>> from mpmath import *
>>> ulp = 2**-54
>>> x = -r + ulp            # -0.36787944117144228
>>> eW(x, g=g0)           # simple version, half precision, as expected.
0 0.367879447562281
1 0.367879447022175
2 0.367879445804813
3 0.367879448020601
4 0.367879446543632
mpf('0.36787944701838976')

>>> eW(x, g=g1)           # log1p version
0 0.367879447562281
1 0.367879447562281
2 0.367879447562281
3 0.367879447562281
4 0.367879447562281
mpf('0.36787944756228136')

>>> eW(x, g=g2)           # log1p + more precise 1/e
0 0.367879447562281
1 0.367879446846838
2 0.367879446801743
3 0.367879446801563
4 0.367879446801563
mpf('0.36787944680156281')

>>> exp(lambertw(x))    # accurate eW
mpf('0.36787944680156281')