(28/48/50) Lambert W Function
|
03-23-2023, 02:56 AM
(This post was last modified: 04-04-2023 01:41 PM by Albert Chan.)
Post: #9
|
|||
|
|||
RE: (28/48/50) Lambert W Function
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) W(-1/e) = -1 --> a ≈ -1/e, x ≈ -1 --> (1+x) is tiny Let 1/e = r + err 1 - ln(x/a) = - ln(x*(r+err)/a) = - log1p((x*err-(a+r)+r*(1+x))/a) x + ln(x/a) = (x+1) - (1 - ln(x/a)) Terms inside log1p sorted in increasing size, r*(1+x) the biggest. With this, and some minor changes to my expW function, I get this: Code: function W(a, x, verbal) lua> W(-0.367, 0, true) -0.9323991847479226 -0.9323991847479282 lua> W(-0.367, -1, true) -1.0707918867680617 -1.0707918867680521 lua> W(1e10, 0, true) 18.111645253927524 20.02372402639066 20.028685384073526 20.02868541330495 20.028685413304952 Update 1: added special case W0(0) = 0, W-1(0) = -∞ (should it be NaN?) Update 2: removed guard "if h == err then return -1 end", see next post. Update 3: improve W guess with a Newton's step, based from y = e^W y = (y+a) / (log(y)+1) → x = a/y = a/(y+a) * (log(y)+1) Update 4: change behavior from "W branch -1 if x" to "x default to branch 0" Assumed input, x = nil/false/0/-1 only: x = nil/false/0 get branch 0; x = -1 get branch -1 |
|||
« Next Oldest | Next Newest »
|
User(s) browsing this thread: 1 Guest(s)