lambertw, all branches - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html) +--- Forum: General Forum (/forum-4.html) +--- Thread: lambertw, all branches (/thread-19784.html) Pages: 1 2 3 lambertw, all branches - Albert Chan - 04-07-2023 01:24 PM I am starting a new thread, to discuss issues of solving W(a, k=0), extended to complex numbers. The thread is based on this post, which supposed guaranteed correct branch, if f=0 (04-02-2023 11:12 PM)Albert Chan Wrote:   (03-26-2023 06:43 PM)Albert Chan Wrote:  Because solving for y = e^W(a) is easier, W code is not recommended. W code is perfect for other branches! (expW code for W branch 0, -1) Numbering the branches of the Lambert W function Except for 2*k*pi*I term, it is basically same formula for W code! lua> f = fn'x: x + I.log(x) - I.log(a) - 2*k*pi*I' lua> df = fn'x: 1 + 1/x' Let X = |x|, A = |a| log(x) - log(a) = log(X/A) + I*(arg(x) - arg(a)) log(X/A) = 1 + log(X*(r+err)/A) ≈ 1 + log1p((X*err - (A+r) + r*(1+X)) / A) This means we can evaluate f accurately around branch point! f=0 guaranteed we get the correct branch, which is nice. There is no need to use complicated W guesses. lua> a, k = 3+4*I, 5 lua> x = 2*k*pi*I -- guess for Wk(a) lua> repeat h=f(x)/df(x); x=x-h; print(x, I.abs(h)) until x == x+1e-6*h (-1.815554248884793+30.714634540472343*I)       1.9462907525577031 (-1.817005891456528+30.71333413897562*I)        0.0019489253984594729 (-1.8170058918466274+30.713334137004896*I)      2.0089622485451688e-009 (-1.8170058918466274+30.713334137004896*I)      0 >>> from mpmath import * >>> lambertw(3+4j, k=5) mpc(real='-1.8170058918466274', imag='30.713334137004896') RE: lambertw, all branches - Albert Chan - 04-07-2023 02:47 PM f = x + ln(x) - ln(a) - 2*k*pi*I Solve for f=0 is easy if x is big. For |k| > 1, we can solve directly, since x imaginery part ≈ 2*k*pi (see OP example) For k=-1, solve for f=0 is especially difficult. Fortunately, we can flip it, to solve for k=1 instead. W(z, k) == conj(W(conj(z), -k) In fact, I flip all negative k's. So now, k is non-negative. (04-02-2023 11:12 PM)Albert Chan Wrote:  W code is perfect for other branches! (expW code for W branch 0, -1) Again, for now on, we assumed k >= 0 f=0 may still have problem converging if k = 0 or 1, especially for small x It does not matter if we have a really good guess. It is why expW code is needed. It handle this issue nicely. (*) expW code also handled accuracy issue elegantly, when a is close to branch point, -1/e Problem is expW code lost phase information, y = e^(x + 2*pi*k*I) = e^x expW can only get W(a, ±1) with the smaller imaginery part. expW code can get this: >>> lambertw(-1-.1j, 1) (-0.373678737258224 + 1.41169809383182j) But not this: >>> lambertw(-1+.1j, 1) (-2.04424331708213 + 7.48779784363809j) Note that negative a imag part corresponded to smaller W(a, k=1) imag part. It can be deduced from f f = x + ln(x) - ln(a) - 2*k*pi*I = 0 x.imag + phase(x) = 2*k*pi + phase(a) z.imag = abs(z) * sin(phase(z)) a.imag < 0      ⇒ phase(a) < 0      ⇒ RHS smaller      ⇒ x.imag smaller (*) Based from tests, solve for y*ln(y)=a is better than x*e^x=a, even with help of Halley's method. mpmath.lambertw make it work only because it throw in extra precisions to compensate. And, it required a lot of code to get a good starting guess. RE: lambertw, all branches - pier4r - 04-07-2023 06:04 PM Side observation. Wouldn't be better to format such posts a bit more and post them in the Article forum? I think they can be more valuable over time if they can get a bit more structured. RE: lambertw, all branches - Albert Chan - 04-07-2023 07:54 PM Hi, pier4r I am still learning, and math is not up to par into article section, where no feedback allowed. Hopefully, someone reading this may suggest a better way. (04-07-2023 02:47 PM)Albert Chan Wrote:  For k=-1, solve for f=0 is especially difficult. Fortunately, we can flip it, to solve for k=1 instead. W(z, k) == conj(W(conj(z), -k) Example, as I am writing code, I realized k "flipping" is unnecessary. With equivalent starting guess, convergence is the same. It is just here, "equivalent" may mean a difference of +0*I vs -0*I Below is work in progress code, without flipping k to non-negative. It simply solve f=0 by Newton's method, with a supplied guess x lua> a, k = -0.1, -1 lua> x = I.log(-a) lua> I.W(a, k, true, x) (-2.3025850929940455+0*I) (-3.7769077194622533-11.106812831380271*I) (-3.949948702089913-0.9503184433887455*I) (-3.558355491208183-0.04104686782028688*I) (-3.577077533853747+8.297706605468053e-005*I) (-3.5770735614781484-8.721290183674583*I) ... lua> I.W(a, k, true, I.conj(x)) -- +0*I to -0*I fixed it! (-2.3025850929940455-0*I) (-3.7769077194622533-0*I) (-3.579124176402325-0*I) (-3.5771522746958055-0*I) (-3.5771520639573-0*I) (-3.577152063957297-0*I) (04-07-2023 02:47 PM)Albert Chan Wrote:  f = x + ln(x) - ln(a) - 2*k*pi*I = 0 x.imag + phase(x) = 2*k*pi + phase(a) Guess x with +0*I is a lousy one, with phase(x) = +pi LHS = 0 + (+pi) = +pi RHS = 2*(-1)*pi + (pi) = -pi guess of conj(x) with -0*I is much better, LHS = 0 + (-pi) = -pi = RHS The fix is simple, argument for log, instead of wrong implicit I.new(-a), use explicit -I.new(a) lua> I.log(-I.new(a)) -- good W-1 guess (-2.3025850929940455-0*I) RE: lambertw, all branches - Albert Chan - 04-07-2023 08:40 PM I am keeping code simple. Edge cases, example W0(0) = 0, are not implemented. Code: ```r, err = 0.36787944117144233, -1.2428753672788363E-17 require'complex' function complex.log1p(x) local y = x-(-1); return I.log(y) - ((y-1-x)/y+0) end function complex.W01(a, K, verbal, y)   -- K = abs(W branch) = 0 or 1     K = K or 0     local h, s = a+r+err, function(y) return I.log1p((y-r-err)/r) end     if not y then         y = I.sqrt(2*r*h)*(1-2*abs(K))  -- assume close to branch point         y = r + y * I.sqrt(1+y/(3*r))   -- assume W imag part tiny     end     for i = 1, 10 do         if verbal then print(a/y, 'W01') end         h = y - (y+a) / s(y)         y = y - h -- Newton's correction         if y == y+h*0x1p-13 or not I.finite(y) then return a/y end     end     return NaN end function complex.W(a, k, verbal, x)     a, k = I.new(a), k or 0     local K, h, s = abs(k), a+r+err     if I.abs(h)<.25 and (K==0 or K==1 and signbit(k*I.imag(a))) then         return complex.W01(a, K, verbal), 'W01'     end     local A, T = I.abs(a), 2*k*pi+I.arg(a)     s = function(x) return x + I.new(log(I.abs(x)/A), I.arg(x)-T) end     if not x then     -- W rough guess         if K > 1 then    x = 2*k*pi*I         elseif K==1 then x = A==1 and 1 or I.log(-a)         else --[[K==0]]  x = A==1 and 1 or I.log1p(a)         end     end     for i = 1, 10 do         if verbal then print(x) end         h = x/(x-(-1)) * s(x)         x = x - h -- Newton's correction         if x == x+h*0x1p-13 or not I.finite(x) then return x end     end     return complex.W01(a, K, verbal, a/I.real(x)), 'W01' end``` W01 is solving W branch 0 or an ambiguous ±1 branch (with the small imag part) The code is solving for y*ln(y) = a --> y = e^W(a) --> W(a) = a/y We can "pull out" conjugate: W01(conj(a),k) = conj(W01(a,k)) lua> I.W(3+4*I, 5, true) (0+31.41592653589793*I) (-1.815554248884793+30.714634540472343*I) (-1.8170058914565277+30.71333413897562*I) (-1.8170058918466272+30.713334137004896*I) (-1.8170058918466274+30.713334137004896*I) lua> I.W(3+4*I, -5, true) (-0-31.41592653589793*I) (-1.7565805546696363-28.861921295458465*I) (-1.7547636407605829-28.85710103800315*I) (-1.7547636279545493-28.85710104741841*I) (-1.7547636279545493-28.85710104741841*I) lua> I.W(-0.3,-1,true) (-1.7824743623491122+0*I)      W01 (-1.7813374873395937+0*I)      W01 (-1.7813370234217047+0*I)      W01 (-1.781337023421628+0*I)        W01 lua> I.W(-0.03,-1,true) (-3.506557897319982-0*I) (-5.261733759227882-0*I) (-5.144793893666586-0*I) (-5.1444827237861634-0*I) (-5.144482721515682-0*I) (-5.144482721515682-0*I) lua> I.W(-.1+1e-100*I, -1, true) (-2.3025850929940455-1e-099*I) (-3.7769077194622533+1.2592576024726166e-099*I) (-3.579124176402325-8.545840728943592*I) (-3.7689826579869887-0.8337217972782689*I) (-3.5535464061313635-0.025256342802877585*I) (-3.5771468383315437+6.535012002211599e-005*I) (-3.5770902410326704-8.721224618630954*I) (-3.7811926340275566-0.8444302692850991*I) (-3.5536175209072614-0.02647513850865879*I) (-3.5771431199160215+6.824715919089008e-005*I) (-3.5770875001635236+0*I)      W01 (-3.5771520641833883-1.3880349425219542e-099*I)      W01 (-3.5771520639572962-1.3880252213229777e-099*I)      W01 (-3.5771520639572976-1.388025221322979e-099*I)        W01 W code f=0 have trouble converging. Switched to W01 help a lot. W01 use W last guess, but dropped the unreliable imaginary part. RE: lambertw, all branches - Albert Chan - 04-08-2023 03:21 PM (04-07-2023 07:54 PM)Albert Chan Wrote:  It is just here, "equivalent" may mean a difference of +0*I vs -0*I Just to show how sensitive this ±0*I business is ... lua> I.W(-1e-9,0,true) (-1.0000000005000001e-009+0*I) (-1.000000001e-009+0*I) (-1.0000000010000002e-009+0*I) lua> I.W(I.conj(-1e-9),0,true) (-1.0000000005000001e-009+0*I) (-1.000000001e-009+6.283185316604365e-009*I) (3.1450880362038825e-008-4.729032266472273e-010*I) (-7.848729544901978e-008-9.717506161329017e-008*I) ... 2nd should be complex conjugate of first. W(-1E-9 ± 0*I) ≈ -1.000000001E-9 ± 0*I Imaginary part sign of (a, W0(a)) should match. At first, I thought it is also a case of bad W guess. To keep x imag part intact, x+1 should be written as x-(-1) I added +0 to correction too, just to be sure. Technically, it is not needed. Signed zero: (±0) - 0 = (±0)      ,       (±0) + 0 = 0 Code: ```< function complex.log1p(x) local y = x+1; return I.log(y) - (y-1-x)/y end > function complex.log1p(x) local y = x-(-1); return I.log(y) - ((y-1-x)/y+0) end``` This get W0 guess phase correct, however iterations still diverges. It is interesting Newton's 1st iteration essentially undo log1p patches. lua> I.W(I.conj(-1e-9),0,true) (-1.0000000005000001e-009-0*I) (-1.000000001e-009+0*I) (-1.0000000010000002e-009+6.283185319745956e-009*I) (3.145088037733907e-008-4.729032324603788e-010*I) ... The problem is Newton correction: 1/slope = x/(x+1). We fix it the same way: Code: ```< h = x/(x+1) * s(x) > h = x/(x-(-1)) * s(x)``` lua> I.W(I.conj(-1e-9),0,true) (-1.0000000005000001e-009-0*I) (-1.000000001e-009-0*I) (-1.0000000010000002e-009-0*I) I had applied above 2 patches to previous post. RE: lambertw, all branches - Albert Chan - 04-08-2023 05:54 PM (04-08-2023 03:21 PM)Albert Chan Wrote:  The problem is Newton correction: 1/slope = x/(x+1). We fix it the same way: Code: ```< h = x/(x+1) * s(x) > h = x/(x-(-1)) * s(x)``` Patch is wrong. It get the correct result by accident. The problem is (±0) + 0 = 0, losing sign information. If h imag part = -0, x = x - h may lose x sign information. We could do x = x - (h + 0), but that is ugly. I am sick of this ±0*I, and simply flip k to avoid it. I.W(I.conj(-1e-9), 0) --> I.conj(I.W(-1e-9, -0)) We flip k to ensure x imag part always non-negative. Losing sign information now becomes an advantage. x = x - h will not generate x with imag part of -0. Only update is complex.W code, the rest remains unaffected. I am leaving 1/slope = x/(x-(-1)), but x/(x+1) will work as well. Code: ```function complex.W(a, k, verbal, x)     a, k = I.new(a), k or 0     local a0, h, s = a, a+r+err     local flip = k<0 or k==0 and signbit(I.imag(a))     flip = flip and I.conj or function(x) return x end     a, k = flip(a), abs(k) --> I.imag(x) non-negative         if I.abs(h)<.25 and (k==0 or k==1 and signbit(I.imag(a))) then         return complex.W01(a0, k, verbal), 'W01'     end     local A, T = I.abs(a), 2*k*pi+I.arg(a)     s = function(x) return x + I.new(log(I.abs(x)/A), I.arg(x)-T) end     if not x then     -- W rough guess         if k > 1 then    x = 2*k*pi*I         elseif k==1 then x = A==1 and 1 or I.log(-a)         else --[[k==0]]  x = A==1 and 1 or I.log1p(a)         end     end     for i = 1, 10 do         if verbal then print(flip(x)) end         h = x/(x-(-1)) * s(x)         x = x - h -- Newton's correction         if x == x+h*0x1p-13 or not I.finite(x) then return flip(x) end     end     return complex.W01(a0, k, verbal, a0/I.real(x)), 'W01' end``` RE: lambertw, all branches - Albert Chan - 04-09-2023 03:59 AM complex.W code version 2, with many improvements. W01 function removed, code placed inside W, for branch point cases only. This allow for brute force confirmation that it always converge to correct branch. Result guaranteed finite, NaN check (to get out of infinite loops) not needed. replaced hexfloat 0x1p-13, 0x1p-26 to more readable 1e-4, 1e-8 W guess improved, for all branches. f formula rearranged, to cause catastrophic cancellation, on purpose! Imag part (if small) are not very good, but at least it converge quickly. We got very good estimated W real part, and a reasonable sized imag part. Added convergence confirmation, based from loop counts. If we reached solution before max loops, we have convergence. If not, result may still be good, but Wk(a) likely purely imaginary.    if |Im(x)| < pi, we are on W branch 0, or branch ±1 with small imag part. To be safe, I lower the limit, from pi to 3. Within this region, we can improve result by y = e^W Newton step. With good W real part, one correction is enough, to fix imag part. y = (a+y) / (ln(y) + 1) a/x = (a+a/x) / (ln(a/x) + 1) x = x/(x+1) * (ln(a/x) + 1)   Note that all iterations are using equivalent f = x - (ln(a) - ln(x)) - 2*k*pi*i W01 code use accurate denominator around branch point, ln(y)+1 = log1p((y-1/e)*e) W code use the general form of f, for all k's, but may cost of some accuracy, if |k|≤1 W final correction (with restriction, see next post) use accurate denominator around 0, ln(y)+1 Code: ```r, err = 0.36787944117144233, -1.2428753672788363E-17 require'complex' function complex.log1p(x) local y = x-(-1); return I.log(y) - ((y-1-x)/y+0) end function complex.W(a, k, verbal)     a, k = I.new(a), k or 0     local flip = function(x) return x end   -- identity function     if k<0 or k==0 and signbit(I.imag(a)) then flip=I.conj; a=flip(a); k=-k end     local h, small_imag = a+r+err, (k==1 and signbit(I.imag(a)))     if I.abs(h)<.25 and (k==0 or small_imag) then         local y = I.sqrt(2*r*h) * (1-2*k)   -- close to branch point         y = r + y*I.sqrt(1+y/(3*r))         -- with small imag part         repeat             if verbal then print(flip(a/y)) end             h = y - (y+a) / I.log1p((y-r-err)/r)             y = y - h -- Newton's correction         until y == y+h*1e-4         return flip(a/y), true     end     local A, T = I.abs(a), (2*k*pi+I.arg(a))*I     local x = T -- W rough guess     if small_imag and A<.6 then x = I.log(-a) end     if k==0 then x = I.abs(a+1)<.6 and T/2 or I.log1p(a) end     k = 40      -- max loops     repeat         if verbal then print(flip(x)) end         h = x/(x+1) * (x + I.new(log(I.abs(x)/A),I.arg(x)) - T)         x, k = x-h, k-1     until k==0 or x == x+h*1e-8 or not I.finite(x)     if I.imag(x) >= 3 then return flip(x), k>0 end     if verbal then print(flip(x)) end     return flip(x/(x+1) * (I.log(a/x)+1)), k>0 end``` lua> I.W(3+4*I, 5, true) (0+32.34322175389954*I) (-1.8166634488635778+30.71625713741827*I) (-1.8170058961963056+30.713334138441088*I) (-1.8170058918466274+30.713334137004896*I)      true lua> I.W(3+4*I, -5, true) (0-30.48863131789632*I) (-1.754507767775116-28.860288696238875*I) (-1.7547636338945982-28.857101048916533*I) (-1.7547636279545493-28.85710104741841*I)      true lua> I.W(-0.3678, 0, true) (-0.9793607152032429+0*I) (-0.9793607149578304+0*I) (-0.9793607149578306+0*I)      true lua> I.W(-.1+1e-100*I,-1,true) (-2.3025850929940455-1e-099*I) (-3.7769077194622533-5.084465616380436e-100*I) (-3.579124176402325-5.180347740277737e-100*I) (-3.5771522746958055-5.181454351722726e-100*I) (-3.5771520639573-5.181454470168089e-100*I) (-3.577152063957297-5.18145447016809e-100*I) (-3.577152063957297-1.3880252213229779e-099*I)      true lua> I.W(-.1+1e-100*I,0,true) (-0.10536051565782631+1.1111111111111111e-100*I) (-0.11161906744205384+1.1848854610650868e-100*I) (-0.11183232962808998+1.1874337731206376e-100*I) (-0.11183255915869776+1.1874365171431562e-100*I) (-0.11183255915896297+1.1874365171463266e-100*I) (-0.11183255915896298+1.259138243719729e-100*I)      true lua> I.W(-.1+1e-100*I,1,true) (0+9.42477796076938*I) (-4.330508133087425+7.394500450320804*I) (-4.448949088733658+7.307107936743056*I) (-4.449098178536345+7.307060789152774*I) (-4.44909817870089+7.3070607892176085*I)      true lua> x = 5*I lua> a = x * I.exp(x) lua> x + x:log() - a:log() (0+6.283185307179586*I) lua> I.W(a, 1) (3.1321400622130957e-068+5*I)      false RE: lambertw, all branches - Albert Chan - 04-09-2023 04:36 PM (04-09-2023 03:59 AM)Albert Chan Wrote:  [*] if |Im(x)| < pi, we are on W branch 0, or branch ±1 with small imag part. (04-07-2023 02:47 PM)Albert Chan Wrote:  f = x + ln(x) - ln(a) - 2*k*pi*I = 0 x.imag + phase(x) = 2*k*pi + phase(a) W code forced k non-negative, to guarantee x.imag ≥ 0 Since x.imag and phase(x) have same sign, we have: k=0: (x.imag < phase(a) < pi)              // for k=0, we forced (0 ≤ phase(a) ≤ pi) k=1: (x.imag < pi) ⇒ (LHS < pi + pi) ⇒ (phase(a) < 0), thus W branch ±1 with small imag part. Quote:To be safe, I lower the limit, from pi to 3. Experiments suggested |Im(x)| < pi is a hard limit. (why?) lua> a = -100-92*I lua> x = I.W(a, 1) lua> x -- abs imag part < pi (3.382739150014567+3.1375382604161515*I) lua> x/(x+1) * (I.log(a/x) + 1) -- y = e^x Newton step (3.382739150014567+3.1375382604161515*I) If |Im(x)| > pi, y = e^W Newton step will not improve, but overshoot toward W branch 0. lua> a = -100-93*I lua> x = I.W(a, 1) lua> x -- abs imag part > pi (3.386390672646167+3.1426530733361844*I) lua> x/(x+1) * (I.log(a/x) + 1) -- y = e^x Newton step (4.064553974086096-2.1939787611722497*I) Hard pi limit puzzle solved! ln(a/x) ≈ ln(e^x) = Re(x) + i * smod(Im(x), 2*pi)         // signed mod to limit arg(z) within ±pi If |Im(x)| < pi, ln(a/x) ≈ ln(e^x) = x → x/(x+1) * (ln(a/x) + 1) ≈ x/(x+1) * (x+1) = x            // Newton's step work as expected If pi < |Im(x)| < 3*pi, ln(a/x) ≈ x ± 2*pi*i                       // sign opposite of Im(x) → x/(x+1) * (ln(a/x) + 1) ≈ x + x/(x+1)*(±2*pi*i)     // Newton's step overshoot toward W branch 0 RE: lambertw, all branches - Albert Chan - 04-10-2023 04:44 PM (04-09-2023 03:59 AM)Albert Chan Wrote:  [*] if |Im(x)| < pi, we are on W branch 0, or branch ±1 with small imag part. Another way to look at this, if above is true, we can simplify f: < f = x - (ln(a) - ln(x)) - 2*k*pi*i > f = x - ln(a/x) f general form has issues with catastrophic cancellation, simplifed f is much better. The problem is we don't know |Im(x)| < pi, until we solve it. Also, f simplifed form have multiple roots. A good guess is needed to get correct branch. For f general form, with 1 root, guess good enough for convergence is all that is needed. We use f general form for *good* guess, correct with using simplified f (if applicable) lua> a = 1e-100-1e-10*I lua> I.W(a, 0, true) (1e-100-1e-010*I) (1.000000082740371e-020-1e-010*I) (1.000000082740371e-020-1e-010*I)        -- general f=0 solution (1.0000000000000002e-020-1e-010*I)      true Here, initial guess log1p(a) ≈ a is not good (I.log code need work ...) But, it doesn't matter. Guess is good enough for general f to converge. f simplified form, with 1 Newton step, we get W(a) ≈ a - a*a RE: lambertw, all branches - Albert Chan - 04-10-2023 06:47 PM (04-10-2023 04:44 PM)Albert Chan Wrote:  Here, initial guess log1p(a) ≈ a is not good (I.log code need work ...) But, it doesn't matter. Guess is good enough for general f to converge. As expected, better complex.log make little difference. I considered Kahan's clog(z), but Python's cmath.log code is simpler. Code: ```function complex.log(z)     local x, y = I.real(z), I.imag(z)     local h = hypot(x, y)     if 0.71 <= h and h <= 1.73 then         if abs(x) < abs(y) then x,y = y,x end         h = log1p((x+1)*(x-1) + y*y)/2     else         h = log(h)     end     return I.new(h, I.arg(z)) end``` lua> a = 1e-100-1e-10*I lua> I.W(a, 0, true) (5.0000000000000005e-021-1e-010*I)      -- log1p(a) ≈ a - a^2/2 (1.0000000413701856e-020-1e-010*I) (1.0000000413701856e-020-1e-010*I) (1.0000000000000001e-020-1e-010*I)      true RE: lambertw, all branches - Albert Chan - 04-13-2023 03:03 PM Here is equivalent complex.W version 2, adapted to machine without signed zero. Without signed zero, this trick is out: W(a, k) = conj(W(conj(a), -k) Instead, let W(a, k) = x = re(z) + s*im(z)*I , such that s = ±1, arg(z) ≥ 0 s = sign(arg(x)) = ±1            |x| = |z|            arg(x) = s * arg(z) lua> I.W(-1e-100, -1, true) (-230.25850929940458-0*I) (-235.72143712473638-0*I) (-235.72115887568603-0*I) (-235.72115887568532-0*I) (-235.72115887568535-0*I)      true Above example, s=-1, all x iterations, arg(x) = s * pi Almost identical, for efficiency: arg(x) ≠ bad = s * -pi In other words, iteration overshoot allowed, where arg(z) ≥ 0 does not hold. But, not to the extreme opposite direction, i.e. arg(z) ≠ -pi We need to add arg(x) edge case, to make W (without signed zero) work the same as I.W Code: ```        h = arg(x)         if h == bad: h = -h``` I use Python mpmath, which does not have signed zero. Code: ```from mpmath import *    # no signed zero support r, err = 0.36787944117144233, -1.2428753672788363E-17 def W(a, k=0, verbal=False):     h, s = a+r+err, im(a)<0     small_imag = k==1 and s or k==-1 and not s     if abs(h)<.25 and (k==0 or small_imag):         y = sqrt(2*r*h) * (1-2*k*k) # close to branch point         y = r + y*sqrt(1+y/(3*r))   # with small imag part         while 1:             if verbal: print a/y             h = y - (y+a) / log1p((y-r-err)/r)             y = y - h             if y == y+h*1e-4: return a/y, True     s = 1 - 2 * (k<0 or k==0 and s)     A, T, bad = abs(a), (2*k*pi+arg(a))*j, -s*pi     x = T   # W rough guess     if small_imag and A<.6: x = log(-a)     if k == 0: x = T/2 if abs(a+1)<.6 else log1p(a)     if s*im(x) < 0: x = conj(x)     # s = sign(im(x))     k = 40  # max loops     while 1:         if verbal: print x         h = arg(x)         if h == bad: h = -h         h = x/(x+1) * (x + mpc(log(abs(x)/A),h) - T)         x, k = x-h, k-1         if k==0 or x==x+h*1e-8 or not isfinite(x): break     if s*im(x) >= 3: return x, k>0     if verbal: print x     return x/(x+1) * (log(a/x)+1), k>0``` >>> W(-1e-100, -1, True) -230.258509299405 (-235.721437124736 + 0.0j) (-235.721158875686 + 0.0j) (-235.721158875685 + 0.0j) (mpc(real='-235.72115887568535', imag='0.0'), True) RE: lambertw, all branches - floppy - 04-13-2023 04:14 PM Nice. Since I see the lambert as usefull for infinite iteration, do you have any idea how to use the lambert (with "jumping" into branches?) for following type of problem a) z = x + i *y b) f(z) = sqrt(x*y) + i* (x+y)/2 c) where finally f(f(f..(z)))) = a + i*a (a=AGM(x,y). AGM arithemic geometric mean) RE: lambertw, all branches - Albert Chan - 04-19-2023 01:30 AM (04-07-2023 02:47 PM)Albert Chan Wrote:  (*) Based from tests, solve for y*ln(y)=a is better than x*e^x=a, even with help of Halley's method. Not when W guess is extremely close to true root. We can use Newton's method, f = x*e^x - a, for finishing touches Since we are now using W loops for guess, we might as well make it slightly worse. Instead of log(abs(x)/A), make it log(abs(x)) - log(A), literally, f = x + ln(x) - lnk(a) Code is easier to read, and cheaper to compute. (last term is a constant) Instead of comparing *both* complex parts converged, we test for relative errors. This solved the problem if any complex parts is bouncing around zero. For symmetry, I used relative error test for Y loop as well. With relative error test, it should not get stuck in infinte loops → loop limits removed complex.W code, version 3 Code: ```r, err = 0.36787944117144233, -1.2428753672788363E-17 require'complex' function complex.W(a, k, verbal)     a, k = I.new(a), k or 0     local flip = function(x) return x end   -- identity function     if k<0 or k==0 and signbit(I.imag(a)) then flip=I.conj; a=flip(a); k=-k end     local h = I(a:real() + r + err, a:imag())     local small_imag = (k==1 and signbit(I.imag(a)))     if I.abs(h)<.25 and (k==0 or small_imag) then         local y = I.sqrt(2*r*h) * (1-2*k)   -- close to branch point         y = r + y*I.sqrt(1+y/(3*r))         -- with small imag part         repeat             if verbal then print(flip(a/y)) end             h = y - (y+a) / I.log1p((y-r-err)/r)             y = y - h         until I.abs(h/y) < 1e-9         return flip(a/y)     end     local A, T = I.abs(a), I.new(0,2*k*pi+I.arg(a))     local x = T     -- W rough guess     if k==0  then x = (A+100)/(A+50); x = I.log1p(a*x) / x end     if small_imag then x = I.log(-a) + (A<.5 and 0 or T/2) end     T = T + log(A)  -- = ln_k(a)     if not finite(A) then return flip(T) end     if A==0 then return flip(k==0 and a or T-pi*I) end     if verbal and verbal~=true then x = flip(verbal) end     repeat         if verbal then print(flip(x)) end         h = x/(x+1) * (x + I.log(x) - T)         x = x - h     until not (I.abs(h/x) >= 1e-6)     if verbal then print(flip(x)) end     h = (x-a*I.exp(-x))/(x+1)   -- final touch     return flip(I.finite(h) and x-h or x) end``` (04-09-2023 03:59 AM)Albert Chan Wrote:  lua> x = 5*I lua> a = x * I.exp(x) lua> x + x:log() - a:log() (0+6.283185307179586*I) lua> I.W(a, 1) (3.1321400622130957e-068+5*I)      false complex.W version 2 was not able to converge after 40 loops. Redo with version 3: lua> a (4.794621373315692+1.4183109273161312*I) lua> I.W(a, 1, true) (0+6.570796326794897*I) (-0.033367069043767406+4.994921913968373*I) (1.97521636071743e-005+5.000010560789135*I) (3.753082249819371e-012+5.000000000009095*I) (1.0141772755203484e-016+5*I) (-6.474631941441244e-017+5*I) We need more precisions to get real part right. For reference, with hexfloat of a: W1(a) ≈ 4.51974e-018 + 5.00000 I Udpate 1/19/24: added |a| = ∞ or 0 edge cases lua> I.W(inf), I.W(-inf) (inf+0*I)      (inf+3.141592653589793*I) lua> z = 0 lua> I.W(-z,-1), I.W(z,-1) (-inf-0*I)      (-inf-3.141592653589793*I) Update 1/29/24 < h = (a + r + err) > h = I(a:real() + r + err, a:imag())      -- keep -0*I sign Update 2/01/24 lyuka branch relative limit reduced, from 1e-12 to 1e-9 Although convergence near -1/e is hard, lyuka (improved) formula is very accurate there. Further away, guess formula is not as good, but convergence O(rel_err^2) constant < 1 I also removed lua implementation of log1p(x). it had moved to C code RE: lambertw, all branches - Albert Chan - 04-23-2023 02:49 PM Previous post, complex.W version 3 old guess code (OK, but now replaced) Note: k flipped to non-negative, to force Im(x) non-negative too, to avoid -0*I issue. Code: ```    local A, T = I.abs(a), (2*k*pi+I.arg(a))*I     local x = T     -- W rough guess     if k==0 then x = I.abs(a+1)<.6 and T/2 or I.log1p(a) end     if small_imag and A<.6 then x = I.log(-a) end``` k=0 and Im(a)≥+0, W0(a) guess of log1p(a) is good, except around log1p singularity, a ≈ -1 If a is real, inside (-1,-1/e), iterations may get stuck on real line, or take a while to get out. That's why I had added a "hole" around -1, with purely complex guess. x + ln(x) = ln(|a|) + T Im(x + ln(x)) = T/i = [pi/2, pi] a ≈ -1      → guess x = T/2 = [pi/4, pi/2]*i We could hide log1p singularity inside branch region, for simpler and better W0 guess Code: ```< if k==0 then x = I.abs(a+1)<.6 and T/2 or I.log1p(a) end > if k==0 then x = I.log1p(2*a)/2 end``` New guess is better for small a. W0(ε)          = ε - ε² + 3/2*ε³ - ... log(1+2ε)/2 = ε - ε² + 4/3*ε³ - ... New guess is worse for huge a. (factor of 2) But, it doesn't matter. x ≫ ln(x), curve is very flat. Cost is at most 1 extra iteration. Update: perhaps we can adjust constant 2, to have best of both worlds. Constant is 2 when a is small, 1 when a is huge (getting back log1p(a)) Code: ```< if k==0 then x = I.log1p(2*a)/2 end > if k==0 then x = (A+100)/(A+50); x = I.log1p(a*x) / x end``` k=1 and Im(a)≤-0, W1(a) guess of log(-a) is good, except a ≈ -1 (log(1)=0, terrible W guess) If a is real < -1/e, iterations may get stuck on real line, or take a while to get out. My solution was only use log(-a) for a ≈ 0, and purely imaginary T as default W guess. T/i = 2*pi + arg(a) re(a)≤0:      T/i = 2*pi .+ (-pi, -pi/2] = (pi, 3/2*pi] re(a)≥0:      T/i = 2*pi .+ [-pi/2, -0.] = [3/2*pi, 2*pi] Another approach is to add reasonable complex part to log(-a). if Re(h = a+1/e) < 0, W is complex. Code: ```< if small_imag and A<.6 then x = I.log(-a) end > if small_imag then x = I.log(-a) + (I.real(h)<0 and T/2 or 0) end``` New guess, Re(a) < -1/e, which covered a ≈ -1: arg(a) = (-pi, -pi/2] im(log(-a) + T/2) = (arg(a)+pi) + (pi+arg(a)/2) = (0, pi/2] + (pi/2, 3/4*pi] = (pi/2, 5/4*pi] I've patched new W guess to complex.W, version 3. Better W guess welcome ... Update: it seems to fit better if we again treat a ≈ 0 as the hole. Code: ```< if small_imag then x = I.log(-a) + (I.real(h)<0 and T/2 or 0) end > if small_imag then x = I.log(-a) + (A<.5 and 0 or T/2) end``` Outside the hole, with Re(a) > 0, we have same formula, but different domain. arg(a) = [-pi/2, -0] im(log(-a) + T/2) = (arg(a)+pi) + (pi+arg(a)/2) = [pi/2, pi] + [3/4*pi, pi] = [5/4*pi, 2*pi] Added T/2 match asymptotes imaginary part better. lua> I.W(1e300-I, 1, true) (690.7755278982137+6.283185307179586*I) -- 2*pi (684.2471217875482+6.274017030794511*I) (684.2471666555265+6.274016338969102*I) (684.2471666555265+6.274016338969102*I) -- ≈ 2*pi lua> I.W(1-1e300*I, 1, true) (690.7755278982137+3.9269908169872414*I) -- 5/4*pi (684.2471409153251+4.705523410032926*I) (684.2471850188219+4.705512170223622*I) (684.2471850188219+4.705512170223621*I) -- ≈ 3/2*pi RE: lambertw, all branches - Albert Chan - 04-23-2023 04:40 PM complex.W version 3, adapted for system without signed zero. This time I made a customized LN(x) to honor signed zero. Code: ```from mpmath import *    # no signed zero support r, err = 0.36787944117144233, -1.2428753672788363E-17 def W(a, k=0, verbal=False):     h, s = a+r+err, im(a)<0     small_imag = k==1 and s or k==-1 and not s     if abs(h)<.25 and (k==0 or small_imag):         y = sqrt(2*r*h) * (1-2*k*k) # close to branch point         y = r + y*sqrt(1+y/(3*r))   # with small imag part         while 1:             if verbal: print a/y             h = y - (y+a) / log1p((y-r-err)/r)             y = y - h             if abs(h/y) < 1e-9: return a/y     s = 1 - 2 * (k<0 or k==0 and s)     def LN(x, bad=-s*pi):         x = ln(x)         return conj(x) if im(x)==bad else x     A, T = abs(a), mpc(0,2*k*pi+arg(a))     x = T          # W rough guess     if k==0:  x = (A+100)/(A+50); x = log1p(a*x) / x     if small_imag: x = LN(-a) + (0 if A<.5 else T/2)     T = T + log(A) # = ln_k(a)     if not isfinite(A): return T     if A==0: return a if k==0 else T - sign(k)*pi*j     while 1:         if verbal: print x         h = x/(x+1) * (x + LN(x) - T)         x = x - h         if not (abs(h/x) >= 1e-6): break     if verbal: print x     h = (x-a*exp(-x))/(x+1) # final touch     return x-h if isfinite(h) else x``` Trivia: W(a < -1/e) is complex, W0(a) = conj(W-1(a)) x + ln(x) = lnk(a) ln0(a) = ln(|a|) - pi*j + 2*pi*j = ln1(conj(a)) W0(a) = W1(conj(a)) = conj(W-1(a)) >>> a = -pi/2 >>> W(a, 0, True) (0.375191803678016 + 1.59508875325509j) (0.0128176527641888 + 1.58759064221342j) (-5.58573510823744e-5 + 1.57084664181356j) (6.00809277877961e-10 + 1.57079632603836j) (5.72587285769658e-17 + 1.5707963267949j) (1.30059124689582e-17 + 1.5707963267949j) >>> W(a, -1, True) (0.451582705289455 - 1.5707963267949j) (0.0200654830289762 - 1.59320451813763j) (-9.57300550001525e-5 - 1.57091376376939j) (1.39642133562001e-9 - 1.5707963231279j) (-3.14444301571107e-17 - 1.5707963267949j) (-5.01152163329281e-17 - 1.5707963267949j) Confirm x * exp(x) = a, both at the same time. (±pi/2*j) * exp(±pi/2*j) = (±pi/2*j) * (±j) = pi/2 * j² = -pi/2 Udpate 1/19/24: added |a| = ∞ or 0 edge cases p2> W(inf), W(-inf) ((+inf + 0.0j), (+inf + 3.14159265358979j)) p2> W(0, -1) (-inf - 3.14159265358979j) Default setup is 0 == +0, that is why we get W-1(+0) = -∞ - pi*I To get W-1(0) = -∞, we need special case for k=-1, 0 == −0 + 0*I Real part -0, Imag part +0 ... it is just too confusing. I am keeping 0 == +0, for consistency. ln(0) = ln(+0) = -∞ --> Wk(+0) = 0 if k==0 else -∞ + (2*k-sign(k))*pi*I Of course, we are free to patch it to make it pick what we wanted. But then, what is Wk(0) in general ? Is 0 still meant +0 + 0*I ? This is a weakness of not having signed zero. Update 2/01/24 lyuka branch relative limit reduced, from 1e-12 to 1e-9 (see post #14) RE: lambertw, all branches - Albert Chan - 01-19-2024 04:14 PM I have added w = Wk(a), for edge cases, |a| = ∞ or 0 (updated to post 14 and 16) For |a| → ∞, we have Wk(a) = lnk(a) see https://www.hpmuseum.org/forum/thread-21172-post-182967.html#pid182967 For |a| → 0, k=0 XCas> series(LambertW(x))      → x - x^2 + 3/2*x^3 - 8/3*x^4 + 125/24*x^5 For |a| → 0, k≠0, we use W branch k definition w + ln(w) = lnk(a) = ln(|a|) + (arg(a) + 2*k*pi)*I Let w = r + θ*I, matching complex parts --> (r + re(ln(r)) ≈ r) = ln(|a|) = -∞ --> (θ + sign(θ)*pi) = (arg(a) + 2*k*pi) --> sign(θ + sign(θ)*pi) = sign(arg(a) + 2*k*pi) --> Since k≠0,   sign(θ) = sign(k) = ±1 Wk(a → 0) = a if k==0 else lnk(a) - sign(k)*pi*I This depends on meaning of arg(a) = atan2(im(a), re(a)) = atan2(±0, ±0) IEEE starndard assumed real part is (ε → 0), but imaginery part is truly 0. IEEE: arg(±0 ± 0*I) → arg(±ε ± 0*I) → ±0 or ±pi Technically, we don't really know arg(0), except its full range (-pi, pi] WA avoid all these, ∞ + (finite number, real or complex) = ∞ WA also don't have signed zero, ±0 ±0*I collapsed to plain 0 WA: Wk(0) = 0 if k==0 else -∞ RE: lambertw, all branches - Albert Chan - 01-20-2024 04:48 PM Newton's method convergence to prove previous post. x * e^x = a               // x = Wk(a) x + ln(x) = lnk(a) = c If c = ±∞ + θ*I, then (|x| → ∞) ⇒ (slope = 1+1/x ≈ 1) Iteration formula is practically equivalent to Newton's method. x = c - ln(x) x1 = c                      // guess x2 = c - ln(x1) x3 = c - ln(x2)          // ... until x's converged if a = ∞*cis(θ), r = re(c) = ∞ c = r + θ*I ≈ |r| * cis(θ/r) ln(c) ≈ ln(|r|) + (θ/r)*I x2 ≈ c - ln(c) ≈ (r-ln(|r|)) + (θ-(θ/r))*I ≈ (r + θ*I) = x1 --> Wk(∞*cis(θ)) = lnk(∞*cis(θ)) if a = 0*cis(θ) and k≠0, r = -∞, sign(θ) = sign(k) c = r + θ*I ≈ |r| * cis(sign(k)*(pi - |θ/r|)) ln(c)   ≈ ln(|r|) + (sign(k)*(pi - |θ/r|))*I  ≈ ln(|r|) + sign(k)*pi*I x2 = c - ln(c) ≈ (r-ln(|r|)) + (θ-sign(k)*pi)*I ≈ r + θ2*I sign(θ2) = sign(θ) = sign(k)      // 1-sided convergence (★) ln(x2) ≈ ln(|r|) + (sign(k)*(pi - |θ2/r|))*I ≈ ln(|r|) + sign(k)*pi*I x3 = c - ln(x2) ≈ c - ln(c) = x2 --> Wk≠0(0*cis(θ)) = lnk(0*cis(θ)) - sign(k)*pi*I (★) Without slope correction, convergence is 1-sided, not crossing discontinuity. lua> c = I(-99999*log(10), -pi)      -- = ln-1(-10^-99999) lua> c (-230256.2067143116-3.141592653589793*I) lua> c - I.log(_) (-230268.55366222185-1.3643899976489848e-05*I) lua> c - I.log(_) (-230268.55371584298-5.925215873503475e-11*I) lua> c - I.log(_) (-230268.5537158432+0*I)          -- W-1(-10^-99999) conjugate Note that final imaginery part should be -0 This is why I.W(a, k) flip to positive k, iterate, then flip back for negative k lua> c = c:conj() lua> c (-230256.2067143116+3.141592653589793*I) lua> c - I.log(_) (-230268.55366222185+1.3643899976489848e-05*I) lua> c - I.log(_) (-230268.55371584298+5.925215873503475e-11*I) lua> c - I.log(_) (-230268.5537158432+0*I) lua> _:conj() (-230268.5537158432-0*I)          -- W-1(-10^-99999) RE: lambertw, all branches - Gil - 01-20-2024 10:52 PM I looked your case k=-1 & k=1, x=-.1 + i ×1E-99. With my initial guess= 'LN(((-.1,1.E-100)+1./e)*(e-\v/2.-1.)+\v/(((-.1,1.E-100)+1./e)*2./e)+1./e)' = (-.11284787341,1.26802807724E-100), the built-in multisolver from HP50 runs endless. I looked in your programs, but could not understand which initial Xo formulae I should use for tiny imaginary part. I tried, from what I guessed from your code, as a guess xo=ln(-a), and it worked. When should I used that guess formula for my initial guess? For example, xo=ln(-a) failed to produce any answer for W0((1.E-100,-.0000000001)): a bad guess. I had to use instead the first formulae to get an answer for W0((1.E-100,-.0000000001)): (9.99905166733E-21,-1.00000000005E-10). And what do you mean by I.log(-a)? Is I.log(-x) a special function or is it i×ln(-x). Could you help me on that matter? Thanks in advance. RE: lambertw, all branches - Albert Chan - 01-21-2024 01:14 AM (01-20-2024 10:52 PM)Gil Wrote:  I looked your case k=-1 & k=1, x=-.1 + i ×1E-99. With my initial guess= 'LN(((-.1,1.E-100)+1./e)*(e-\v/2.-1.)+\v/(((-.1,1.E-100)+1./e)*2./e)+1./e)' = (-.11284787341,1.26802807724E-100), the built-in multisolver from HP50 runs endless. x = Wk(a) Newton's method required a good guess to start. (iterations cannot cross discontinuity!) lyuka's e^W formula is only used when a ≈ -1/e, and (k=0 or k=±1 and x with small imag part) Above example probably too far away from -1/e (my code use lyuka if |a+1/e| < 0.25) Also, your formula is only good for k=0. For k=-1, you need to flip sqrt sign. lua> a, e = I(-0.1, 1e-100), exp(1) lua> I.log((a+1/e)*(e-sqrt(2)-1) + I.sqrt((2/e)*(a+1/e)) + 1/e) -- x guess for k=0 (-0.1128478734106899+1.268028077236754e-99*I) lua> I.log((a+1/e)*(e-sqrt(2)-1) − I.sqrt((2/e)*(a+1/e)) + 1/e) -- x guess for k=-1 (-5.225138594746179-9.75118016335884e-98*I) If x imag part is big, it just use that. -- local A, T = I.abs(a), I.new(0,2*k*pi+I.arg(a)) -- local x = T -- W rough guess lua> I.W(a, 1, true) (0+9.42477796076938*I) (-4.330508133087424+7.394500450320804*I) (-4.448949088733658+7.307107936743056*I) (-4.449098178536345+7.307060789152774*I) (-4.44909817870089+7.3070607892176085*I) (-4.44909817870089+7.3070607892176085*I) If k=-1 and x imag part is small, we wanted real guess when a is negative real. Also, x→-∞ if a→0, log(-a) fit the bill. But if a is too negative (left of -1/e), code added some imaginery part. Imag part is needed because a=-1 --> log(-a) = log(1) = 0, very bad x guess! -- if small_imag then x = I.log(-a) + (A<.5 and 0 or T/2) end lua> I.W(a, -1, true) (-2.3025850929940455-1e-99*I) (-3.776907719462253-5.084465616380437e-100*I) (-3.5791241764023245-5.180347740277738e-100*I) (-3.577152274695805-5.181454351722727e-100*I) (-3.5771520639572993-5.18145447016809e-100*I) (-3.577152063957297-1.388025221322979e-99*I)