(28/48/50) Lambert W Function

[Albert Chan]

(03-27-2023 08:57 PM)Albert Chan Wrote:  lua> h = a + r + err
lua> x = sqrt(2*h/r) -- rough guess for e*h = (1+x)*log1p(x) - x
lua> x = x * sqrt(1+x/3) -- better guess
lua> y = r + r*x -- = e^W(a), if x is correct

Technically, it should be y = (1+x)/e = (1+x)*(r+err)
Numerically, err < 1/2 ULP of r, thus y = r + r*x

e*h = (e*y) * log(e*y) - (e*y-1)
h = y * (log(y)+1) - y + 1/e
h - 1/e = y*log(y) = a

---

A novel way is to directly solve for x, using log1p_sub(x), then convert to W(a) or e^W(a)
Except for getting accurate h, err = 1/e - float(1/e) is not used anymore.

Let H=e*h, L = log(1+x) - x, accurately computed

f = (1+x)*log(1+x) - x - H
f' = log(1+x) + (1+x)/(1+x) - 1 = log(1+x)

x - f/f' = (H-L) / (x+L)

Note that numerator is really an addition, with both H, -L positive.

lua> x -- better guess, from above quote
0.020853747742736108
lua> H = h/r
lua> L = log1p_sub(x)
lua> x = (H-L) / (x+L) -- newton step
lua> x -- converged to full precision
0.020853747998545925
lua> r+r*x, log1p(x)-1 -- = e^W(a), W(a)
0.3755511063314775      -0.9793607149578305