(28/48/50) Lambert W Function

[Albert Chan]

Hi, Gil

FNL(x) = accurate ln(1+x)-x, was designed for real number.
It may need work to extend for complex numbers.

Also, it is defined recursively, until x is small, not a IF THEN ELSE structure.
see thread: Accurate x - log(1+x) (for consisteny, I switched order, and do log1p_sub(x))

Anyway, accurate FNL(x) is not needed here.
Here is a direct replacement snippet for my Python W code

I am using formula to get log1p, instead of calling built-in log1p

Quote:y=1+x --> log1p(x) ≈ ln(y) - (y-1-x)/y

Because of inaccurate x-log1p(x), x does not converge, but (1+x) does.
(again, you may forget about checking convergence, and just loop a few times)

Code:

    if abs(h)<.25 and (k==0 or small_imag):
        h = h/r
        x = sqrt(2*h) * (1-2*k*k)   # close to branch point
        x *= sqrt(1 + x/3)
        while 1:
            if verbal: print a/(r+r*x)
            y = 1+x
            z = ln(y) - (y-1-x)/y   # = log1p(x)
            x = (x-z+h)/z
            z = 1+x
            if abs((z-y)/z) < 1e-12: return a/(r+r*x)

With identical termination condition, results compared well to W original version.

Code:

    if abs(h)<.25 and (k==0 or small_imag):
        y = sqrt(2*r*h) * (1-2*k*k) # close to branch point
        y = r + y*sqrt(1+y/(3*r))   # with small imag part
        while 1:
            if verbal: print a/y
            h = y - (y+a) / log1p((y-r-err)/r)
            y = y - h
            if abs(h/y) < 1e-12: return a/y