lambertw, all branches

[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