lambertw, all branches

[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.