lambertw, all branches

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