lambertw, all branches

[Albert Chan]

Hi, Gil

If you want to implement W to HP50g, ignore lua I.W code, and use only Python code (post 16)
My guess is HP50g does not support signed zero anyway.

Python code, including blank lines, is only 34 lines long.
It is all in the code, but this is the gist of it ...

If |k| > 1, just use imaginery part as guess. There is no Newton's over-shoot issue.
If |k| = 1, and k and im(a) has same sign, again, use T as guess

A, T = abs(a), mpc(0,2*k*pi+arg(a))
x = T            # W rough guess

Now the hard part, we need good guess to avoid Newton over-shoot crossed discontinuity.

If |a+1/e| < 0.25, use lyuka's formula and y = e^W iteration formula, return a/y = x
This take cared of singularity at x = -1/e, W₀(-1/e) = W_-1(-1/e) = -1

If |k| = 1, guess x = LN(-a) + (0 if A<.5 else T/2)

We need customized LN() because we lost signed zero information.
If a ≈ -1, LN(-a) is very tiny, possibly 0, thus the need for adding T/2

If k = 0, guess x = log1p(a*s) / s, where s = (A+100)/(A+50) = (2 if A=0, 1 if A=∞)
We have good W0 asymptote estimate when A=0, and A=∞

repeat until loop does Newton's iteration. x = x - f/f'
f is designed to cause catastrophic cancellation, making x possibly not accurate.

Final touch with g = x*e^x - a, g' = (x+1)*e^x
--> h = g/g' = (x - a*e^-x) / (x+1)

Because of the exponential, h may not be finite.
If h is finite, we update for better x; if not, last x is the best we've got.

Quote:You write
Wk(+0) = 0 if k==0 else -∞ + (2*k-sign(k))*pi*I
OK

You write further
Wk(a → 0) = a if k==0 else lnk(a) - sign(k)*pi*I

Wk(+0) is referring Python code, which does not support signed zero.

ln_k(a) = ln(a) + 2*k*pi*I = ln(|a|) + (arg(a) + 2*k*pi)*I

ln_k(0) = ln(|0|) + (arg(0) + 2*k*pi)*I = -∞ + 2*k*pi*I
ln_k(0) - sign(k)*pi*I = -∞ + (2*k-sign(k))*pi*I

With system that support signed zero, 0 may mean (±0 + ±0*I)
This is why ln_k(a) is left unevaluated. arg(0) may have multiple values.

lua> z = 0
lua> I(z,z):arg(), I(z,-z):arg()
0                                  -0
lua> I(-z,z):arg(), I(-z,-z):arg()
3.141592653589793      -3.141592653589793