lambertw, all branches

[Albert Chan]

(04-08-2023 03:21 PM)Albert Chan Wrote:  The problem is Newton correction: 1/slope = x/(x+1). We fix it the same way:

Code:

< h = x/(x+1) * s(x)
> h = x/(x-(-1)) * s(x)

Patch is wrong. It get the correct result by accident.
The problem is (±0) + 0 = 0, losing sign information.
If h imag part = -0, x = x - h may lose x sign information.

We could do x = x - (h + 0), but that is ugly.
I am sick of this ±0*I, and simply flip k to avoid it.

I.W(I.conj(-1e-9), 0) --> I.conj(I.W(-1e-9, -0))

We flip k to ensure x imag part always non-negative.
Losing sign information now becomes an advantage.
x = x - h will not generate x with imag part of -0.

Only update is complex.W code, the rest remains unaffected.
I am leaving 1/slope = x/(x-(-1)), but x/(x+1) will work as well.

Code:

function complex.W(a, k, verbal, x)
    a, k = I.new(a), k or 0
    local a0, h, s = a, a+r+err
    local flip = k<0 or k==0 and signbit(I.imag(a))
    flip = flip and I.conj or function(x) return x end
    a, k = flip(a), abs(k) --> I.imag(x) non-negative    
    if I.abs(h)<.25 and (k==0 or k==1 and signbit(I.imag(a))) then
        return complex.W01(a0, k, verbal), 'W01'
    end
    local A, T = I.abs(a), 2*k*pi+I.arg(a)
    s = function(x) return x + I.new(log(I.abs(x)/A), I.arg(x)-T) end
    if not x then     -- W rough guess
        if k > 1 then    x = 2*k*pi*I
        elseif k==1 then x = A==1 and 1 or I.log(-a)
        else --[[k==0]]  x = A==1 and 1 or I.log1p(a)
        end
    end
    for i = 1, 10 do
        if verbal then print(flip(x)) end
        h = x/(x-(-1)) * s(x)
        x = x - h -- Newton's correction
        if x == x+h*0x1p-13 or not I.finite(x) then return flip(x) end
    end
    return complex.W01(a0, k, verbal, a0/I.real(x)), 'W01'
end