Post Reply 
(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
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
Find all posts by this user
Quote this message in a reply
Post Reply 


Messages In This Thread
(28/48/50) Lambert W Function - John Keith - 03-20-2023, 08:43 PM
RE: (28/48/50) Lambert W Function - Albert Chan - 03-27-2023 11:30 PM
RE: (28/48/50) Lambert W Function - Gil - 01-29-2024, 11:04 AM
RE: (28/48/50) Lambert W Function - Gil - 01-29-2024, 02:47 PM
RE: (28/48/50) Lambert W Function - Gil - 01-29-2024, 06:46 PM
RE: (28/48/50) Lambert W Function - Gil - 01-29-2024, 09:50 PM
RE: (28/48/50) Lambert W Function - Gil - 01-30-2024, 12:33 AM
RE: (28/48/50) Lambert W Function - Gil - 01-30-2024, 12:04 PM
RE: (28/48/50) Lambert W Function - Gil - 01-30-2024, 02:52 PM
RE: (28/48/50) Lambert W Function - Gil - 01-31-2024, 07:10 PM



User(s) browsing this thread: 1 Guest(s)