(28/48/50) Lambert W Function
|
03-26-2023, 06:43 PM
Post: #12
|
|||
|
|||
RE: (28/48/50) Lambert W Function
(03-23-2023 02:56 AM)Albert Chan Wrote: Here is my implementation of accurate x = W(a), both branches. 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 |
|||
« Next Oldest | Next Newest »
|
User(s) browsing this thread: 2 Guest(s)