(28/48/50) Lambert W Function
|
03-27-2023, 11:30 PM
(This post was last modified: 03-30-2023 03:58 PM by Albert Chan.)
Post: #15
|
|||
|
|||
RE: (28/48/50) Lambert W Function
(03-27-2023 08:57 PM)Albert Chan Wrote: lua> h = a + r + err 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 |
|||
« Next Oldest | Next Newest »
|
User(s) browsing this thread: 1 Guest(s)