(28/48/50) Lambert W Function

[Albert Chan]

(03-23-2023 02:56 AM)Albert Chan Wrote:  Here is my implementation of accurate x = W(a), both branches.

x * e^x = a
ln(x) + x = ln(a)      → f = x + ln(x/a) = 0

x - f/f' = x - (x + ln(x/a)) * x/(1+x) = (1 - ln(x/a)) * x/(1+x)

I just realized this is equivalent to solving for y = e^W(a)
With "same" guess, x*y=a, Newton convergence rate are identical.

x = (1 - ln(x/a)) * x/(1+x)
a/y = (1 + ln(y)) * 1/(1+y/a)

y = (y+a) / (1 + ln(y))

lua> expW(1e99, nil, true)
2.5e+098      4
5.492824331831047e+096      182.05570387623288
4.493790313227532e+096      222.52929716290643
4.493356750520384e+096      222.55076895111614
4.493356750426822e+096      222.55076895575016
4.493356750426821e+096      222.55076895575021

lua> W(1e99, nil, true)
182.05570387623249
222.52929716290646
222.55076895111617
222.55076895575016
222.5507689557502

Note: W guess used 1 Newton step. Both code converged the same way.
Because solving for y = e^W(a) is easier, W code is not recommended.

Bonus: with y, there are 2 ways to recover W(a) = a/y or log(y)
Depends on size of y (use log if y is big), we can make W(a) estimate better.

Example, above last 2 y's differ by 1 ULP, but the log's are the same.
(Not exactly. Here, y 1 ULP error translated to 0.007 ULP error in x)

lua> log(4.493356750426822e+096)
222.5507689557502
lua> log(4.493356750426821e+096)
222.5507689557502