lambertw, all branches

[Albert Chan]

(04-09-2023 03:59 AM)Albert Chan Wrote:  [*] if |Im(x)| < pi, we are on W branch 0, or branch ±1 with small imag part.

(04-07-2023 02:47 PM)Albert Chan Wrote:  f = x + ln(x) - ln(a) - 2*k*pi*I = 0
x.imag + phase(x) = 2*k*pi + phase(a)

W code forced k non-negative, to guarantee x.imag ≥ 0
Since x.imag and phase(x) have same sign, we have:

k=0: (x.imag < phase(a) < pi)              // for k=0, we forced (0 ≤ phase(a) ≤ pi)
k=1: (x.imag < pi) ⇒ (LHS < pi + pi) ⇒ (phase(a) < 0), thus W branch ±1 with small imag part.

Quote:To be safe, I lower the limit, from pi to 3.

Experiments suggested |Im(x)| < pi is a hard limit. (why?)

lua> a = -100-92*I
lua> x = I.W(a, 1)
lua> x -- abs imag part < pi
(3.382739150014567+3.1375382604161515*I)
lua> x/(x+1) * (I.log(a/x) + 1) -- y = e^x Newton step
(3.382739150014567+3.1375382604161515*I)

If |Im(x)| > pi, y = e^W Newton step will not improve, but overshoot toward W branch 0.

lua> a = -100-93*I
lua> x = I.W(a, 1)
lua> x -- abs imag part > pi
(3.386390672646167+3.1426530733361844*I)
lua> x/(x+1) * (I.log(a/x) + 1) -- y = e^x Newton step
(4.064553974086096-2.1939787611722497*I)

---

Hard pi limit puzzle solved!

ln(a/x) ≈ ln(e^x) = Re(x) + i * smod(Im(x), 2*pi)         // signed mod to limit arg(z) within ±pi

If |Im(x)| < pi, ln(a/x) ≈ ln(e^x) = x
→ x/(x+1) * (ln(a/x) + 1) ≈ x/(x+1) * (x+1) = x            // Newton's step work as expected

If pi < |Im(x)| < 3*pi, ln(a/x) ≈ x ± 2*pi*i                       // sign opposite of Im(x)
→ x/(x+1) * (ln(a/x) + 1) ≈ x + x/(x+1)*(±2*pi*i)     // Newton's step overshoot toward W branch 0