lambertw, all branches

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