[HP41] Lambert function RPN; question
12-30-2022, 11:48 PM (This post was last modified: 04-04-2023 01:37 PM by Albert Chan.)
Post: #16
 Albert Chan Senior Member Posts: 2,686 Joined: Jul 2018
RE: [HP41] Lambert function RPN; question
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
 « Next Oldest | Next Newest »

 Messages In This Thread [HP41] Lambert function RPN; question - floppy - 12-26-2022, 12:13 PM RE: [HP41] Lambert function RPN; question - Didier Lachieze - 12-26-2022, 02:07 PM RE: [HP41] Lambert function RPN; question - floppy - 12-26-2022, 05:46 PM RE: [HP41] Lambert function RPN; question - Albert Chan - 12-26-2022, 04:04 PM RE: [HP41] Lambert function RPN; question - Sylvain Cote - 12-26-2022, 04:32 PM RE: [HP41] Lambert function RPN; question - floppy - 12-26-2022, 06:45 PM RE: [HP41] Lambert function RPN; question - Thomas Klemm - 12-27-2022, 02:31 AM RE: [HP41] Lambert function RPN; question - Albert Chan - 12-27-2022, 04:26 AM RE: [HP41] Lambert function RPN; question - Thomas Klemm - 12-27-2022, 04:36 PM RE: [HP41] Lambert function RPN; question - Thomas Klemm - 12-27-2022, 04:48 PM RE: [HP41] Lambert function RPN; question - floppy - 12-27-2022, 06:42 PM RE: [HP41] Lambert function RPN; question - Thomas Klemm - 12-27-2022, 05:06 PM RE: [HP41] Lambert function RPN; question - Pekis - 12-28-2022, 11:14 AM RE: [HP41] Lambert function RPN; question - Albert Chan - 12-28-2022, 10:26 PM RE: [HP41] Lambert function RPN; question - Albert Chan - 12-30-2022, 11:27 PM RE: [HP41] Lambert function RPN; question - Albert Chan - 12-30-2022 11:48 PM

User(s) browsing this thread: 1 Guest(s)