lambertw, all branches

[Albert Chan]

Note that W curve is only in 1st and 3rd quadrant.
If x is in ℝ (implied a in ℝ too), they have same sign, then ln(x/a) = ln(|x|) - ln(|a|)

This give me an idea!

Code:

def LN(x, bad=-s*pi):
    x = ln(x)
    return conj(x) if im(x)==bad else x

My Python code used customized LN, to track signed zero, to give correct imag parts.
Example, to simulate ln(-1-0*I), set s=-1 --> LN(-1) = conj(pi*j) = -pi*j

But, what if re(a) were never negative to begin with?
Then, LN is not needed! Unlike ±pi, ±0 are numerically the same thing!

Let s = -1 if re(a)<0 else 1

x * e^x = a
(s*x) * e^x = s*a
ln(s*x) + x = ln(s*a) = B

p2> f = lambda x: show(x) + ln(s*x) - B
p2> df = lambda x: 1 + 1/x
p2> def show(x): print nstr(x, n=10); return x
...
p2> a = mpc(-0.1)
p2> s = -1 if re(a)<0 else 1
p2> B = ln(s*a)

p2> W(a, -1)
(-3.5771520639573 + 0.0j)
p2> findroot(f, x0=B, df=df, solver='newton')
(-2.302585093 + 0.0j)
(-3.776907719 + 0.0j)
(-3.579124176 + 0.0j)
(-3.577152275 + 0.0j)
(-3.577152064 + 0.0j)
(-3.577152064 + 0.0j)
(-3.5771520639573 + 0.0j)

This is nothing special, iterations are exactly the same as W code.
But let's try a *bad* guess, x0=B+j, crossed discontinuity!
(this guess for for W(a,-1) would get stuck in endless loops.)

p2> findroot(f, x0=B+j, df=df, solver='newton')
(-2.302585093 + 1.0j)
(-3.44870968 - 0.2167163462j)
(-3.575061443 + 0.002993043143j)
(-3.577151814 - 6.790747928e-7j)
(-3.577152064 + 1.841873128e-14j)
(-3.577152064 - 4.31950186e-29j)
(-3.577152064 + 0.0j)
(-3.5771520639573 + 0.0j)

We can use better solver, without worrying about over-shooting!

p2> findroot(f, x0=B, df=df, solver='halley')
(-2.302585093 + 0.0j)
(-3.48604165 + 0.0j)
(-3.577141461 + 0.0j)
(-3.577152064 + 0.0j)
(-3.577152064 + 0.0j)
(-3.5771520639573 + 0.0j)

p2> findroot(f, x0=B+j, df=df, solver='halley')
(-2.302585093 + 1.0j)
(-3.730387658 - 0.002761341752j)
(-3.577193828 - 2.173902656e-6j)
(-3.577152064 - 1.483011857e-16j)
(-3.577152064 - 2.350988702e-38j)
(-3.577152064 + 0.0j)
(-3.5771520639573 + 0.0j)

Comment:

What we really wanted is s = -1 if re(x)<0 else 1
This keep |arg(s*x)| closer to 0 then pi.

It happens that if final x is real, sign(x) = sign(a)
Unfortunately, for complex a, we don't know what s is ...