lambertw, all branches

[Albert Chan]

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 W_k(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')

[Albert Chan]

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.

[pier4r]

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.

[Albert Chan]

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)

[Albert Chan]

I am keeping code simple. Edge cases, example W₀(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.

[Albert Chan]

(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, W₀(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 W₀ 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.

[Albert Chan]

(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

[Albert Chan]

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 W_k(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

[Albert Chan]

(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

[Albert Chan]

(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

[Albert Chan]

(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

[Albert Chan]

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)

[floppy]

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)

[Albert Chan]

(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) - ln_k(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: W₁(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

[Albert Chan]

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

[Albert Chan]

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, W₀(a) = conj(W_-1(a))

x + ln(x) = ln_k(a)

ln₀(a) = ln(|a|) - pi*j + 2*pi*j = ln₁(conj(a))

W₀(a) = W₁(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) = -∞

--> W_k(+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 W_k(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)

[Albert Chan]

I have added w = W_k(a), for edge cases, |a| = ∞ or 0 (updated to post 14 and 16)

For |a| → ∞, we have W_k(a) = ln_k(a)
see https://www.hpmuseum.org/forum/thread-21...#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) = ln_k(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

W_k(a → 0) = a if k==0 else ln_k(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: W_k(0) = 0 if k==0 else -∞

[Albert Chan]

Newton's method convergence to prove previous post.

x * e^x = a               // x = W_k(a)
x + ln(x) = ln_k(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

--> W_k(∞*cis(θ)) = ln_k(∞*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

--> W_k≠0(0*cis(θ)) = ln_k(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)

[Gil]

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.

[Albert Chan]

(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 = W_k(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)