[HP41] Lambert function RPN; question

[Albert Chan]

The reason for accurate eW estimate at a ≈ -1/e, is that it take a lot of iterations for convergence.

Once we get close, convergence speed up dramatically (slope = ln(1+Z) ≈ Z)
For IEEE double, if we get around 40 bits, Newton's step converged in 1 step.

a ≈ -1/e --> y ≈ 1/e, and both are inputs (considered exact) --> (y+a) is exact.

We do need to use log1p() instead of ln() for accurate slope.

ln(y)+1 = ln(e*y) = log1p(e*(y-1/e))

Term (1/e) need extra precisions, to avoid cancellation errors.
--> (r + err) = rounded 108 bits of 1/e (about 33 decimals digits accuracy)

Code:

function expW(a, y, verbal)
    y = y or 0  -- default branch 0
    local h, s
    if a <= -0.1 then
        h = a + r + err
        y = sqrt(2*r*h) * (y+y+1)
        y = y * sqrt(1+y/(3*r)) + r
        s = function(y) return log1p((y-r-err)/r) end
    else
        y = (y+1) + (y+y+1)*a/4
        s = function(y) return log(y) + 1 end
    end
    repeat
        if verbal then print(y, a/y) end
        h = y - (y+a) / s(y)
        y = y - h   -- Newton's correction
    until y == y+h*0x1p-13 or y ~= y
    return y, a/y   -- e^W(a), W(a)
end

Convergence should take at most 6 iterations
Accuracy should be +/- a few ulps, for a ≥ -1/e

lua> expW(1e300) -- default = branch 0
1.4614601088436296e+297      684.2472086297608
lua> _ * log(_)
1e+300

lua> expW(-1e-300, -1) -- branch -1
1.4340561272249246e-303      -697.3227762954601
lua> _ * log(_)
-1e-300

Update 1: I implemented Lua code for x = W(a), solving for f = x + ln(x/a) = 0
This post is updated to match that, with a verbal option.
Also, guard removed, "if h == err then return r, -1 end"

Update 2: change behavior from "W branch -1 if y" to "y default to branch 0"
Assumed input, y = nil/false/0/-1 only: y = nil/false/0 get branch 0; y = -1 get branch -1

lua> expW(-0.367, 0, true)
0.3936082265591097              -0.9323992112875368
0.3936082377626915              -0.9323991847479225
0.3936082377626891              -0.9323991847479282