(28/48/50) Lambert W Function

[John Keith]

NOTE: The first and last programs have been obsoleted by the program in post #22, which should be used in most cases. That program gives accurate results over the entire complex plane for any branch, positive or negative, but accuracy is limited for values very close to the branch point -1/e.

For very accurate results near the branch point in branches 0 and -1, use the program in post #28, which is a translation of Albert Chan's HP-71 program in post # 19.

The following text refers to the (mostly obsolete) programs in this post only.

---------------------------------------------------------------

This program computes the Lambert W function W(z) for branches 0 and -1. For branch -1, the program covers the real-valued range -1/e < z < 0. For branch 0, the program covers the entire complex plane. For most values, the results are accurate to within 2 ULPs. For values very close to 1/e, the results are approximate. The program expects the branch (0 or -1) on level 2 and z on level 1.

Some notes and background information:

The gold standard for Lambert W programs on the HP 50g is Egan Ford's SPECIAL.LIB. It is very accurate for almost any input value, and can return results for any branch. It does have some problems, however. It is an HPGCC program and only works on the 49g+ and 50g. It requires two large programs on the calculators SD card. Also, it is rather slow. I believe that the slowness is caused by the need to fetch the executable from the SD card, as one would expect the HPGCC code to be quite fast.

The program described here is comparatively limited since no branches other than 0 and -1 are covered, and then only the real-valued range of branch -1. However, it is about 5 times as fast as SPECIAL.LIB and about 9 times as fast as programs using the calculator's built-in solvers. It also runs on all RPL calculators.

If anyone here has an idea of how to get an accurate estimate for other branches, or to improve accuracy around the branch point at -1/e, please chime in!

The code that provides the initial estimate for branch 0 comes from page 8 of Iacono and Boyd. The estimation code for branch -1 comes from formula 27 on page 12 of https://www.sciencedirect.com/science/ar...via%3Dihub. The recurrence used in the main loop occurs in both papers and is mentioned in Wikipedia. Many additional references here.

Code:


\<< 1. EXP ROT ROT SWAP
  IF NOT                                             @ Branch 0
  THEN DUP ABS .00000001
    IF \>=
    THEN DUP ROT * 1. + \v/ DUP 1.14956131 * 1. +    @ Estimate for
      SWAP 1. + LN .4549574 * 1. + / LN 2.036 * 1. - @ branch 0, |z| > 1E-8
    ELSE DUP ROT DROP                                @ DUP If |z| < 1E-8
    END
  ELSE                                               @ Branch -1
    IF DUP -.25 \<=                                  @ Estimate for
    THEN DUP -1. SWAP 4. ROLL * 1. + 2. * \v/ -      @ z <= 0.25
    ELSE DUP ROT DROP NEG LN DUP NEG LN -            @ For z > 0.25
    END
  END 1. 3.                                          @ Main loop
  START DUP2 / LN 1. + SWAP DUP 1. + / *
  NEXT SWAP DROP
\>>

Also, a streamlined version for the 49 and 50, 25 bytes smaller:

Code:


\<< 1. EXP UNROT SWAP
  IF NOT
  THEN DUP ABS .00000001 \>=
    { DUP ROT * 1. + \v/ DUP 1.14956131 * 1. +
      SWAP 1. + LN .4549574 * 1. + / LN 2.036 * 1. - }
    { NIP DUP } IFTE
  ELSE DUP -.25 \<=
    { DUP -1. SWAP 4. ROLL * 1. + 2. * \v/ - }
    { NIP DUP NEG LN DUP NEG LN - } IFTE
  END 1. 3.
  START DUP2 / LN 1. + SWAP DUP 1. + / *
  NEXT NIP
\>>

Now, here are a couple of experimental programs for those interested in further exploration. The first one is a variation of the recurrence from the main program using the LongFloat Library. The program expects z on level 2 and an estimate (such as the result of the main program) on level 1. It is rather slow, about 22 seconds on my 50g but reasonable on emulators. The value of c assumes the default value of DIGITS which is 50. Removing the R\<-\->F at the end will leave a LongFloat number (type 27) on the stack.

Code:


\<< R\<-\->F SWAP R\<-\->F SWAP
    1. R\<-\->F 1.E-24 R\<-\->F \-> f1 c
  \<<
    DO DUP2 FDIV FLN f1 FADD OVER DUP f1 FADD FDIV FMULT
    UNTIL DUP ROT FSUB FABS c F.LE
    END NIP R\<-\->F
  \>>
\>>

The last program implements the Halley recurrence shown at the Wikipedia link and in the Lóczi paper. The reason for including it is that the recurrence used in the main program works poorly for branches other than branch 0 and the real-valued part of branch -1. It tends to "drag" the values back to branch 0 even if given an accurate estimate for the value in another branch. The Halley recurrence converges rapidly even if the estimate is not very close, although it can be fooled if the estimate is vague. As above, the program expects z on level 2 and an estimate (user-supplied) on level 1. For HP-28 and 48, change PICK3 to 3 PICK.

Code:


\<<
  WHILE DUP2 DUP EXP
    DUP ROT *
    ROT -
    PICK3 1. + ROT *
    DUP2 / ABS .00000000005 > @ 5E-11
  REPEAT PICK3 2. + PICK3 *
    4. PICK 2. * 2. + /
    - / -
  END DROP2
\>>

[Albert Chan]

The issue is how to accurately calculate Newton's iteration, around branch point, -1/e
I prefer solving for y = e^W(a). Its Newton's method is very simple.

x*e^x = ln(y)*y = a
f = y*ln(y) - a

y = y - f/f' = y - (y*ln(y) - a) / (ln(y)+1)
y = (y+a) / (ln(y)+1)

Around a ≈ -1/e, (y, a) are considered input, thus "exact".
They are about the same size, with opposite sign.
By Sterbenz lemma, (y+a) is exact.

e^W(-1/e) = 1/e      // a ≈ -1/e --> y ≈ +1/e

(12-30-2022 11:48 PM)Albert Chan Wrote:  ln(y)+1 = ln(e*y) = log1p(e*(y-1/e))

For accurate denominator, we need "doubled" precision for 1/e
For 12 digits, 1/e ≈ 0.367879441171 + 0.442321595524E-12

To reduce iterations, we also need good starting estimate.
Lua e^W code (0, -1 branch): https://www.hpmuseum.org/forum/thread-19...#pid167919

[Gerald H]

This programme works for positive reals - I'm interested in the comparison of accuracy:

Code:

::
  CK1&Dispatch
  BINT1
  ::
    '
    ID 0
    TRUE
    3PICK3PICK
    DUP'
    %EXP
    '
    %*
    '
    %-
    BINT6
    SYMBN
    3PICK
    5PICK
    DUP
    % 500.
    %<=
    ITE
    ::
      %1+
      %LN
      DUP
      % .0195
      %*
      %1+
      %*
      % .665
      %*
      % .04
      %+
    ;
    ::
      DUP
      %4
      %-
      %LN
      SWAP
      %LN
      %1
      OVER
      %1/
      %-
      SWAP
      %LN
      %*
      %-
    ;
    NUMSOLVE
    2DROP
    SWAP
    ?PURGE_HERE
    SWAP
    xDROP
  ;
;

[John Keith]

(03-20-2023 10:51 PM)Albert Chan Wrote:  For accurate denominator, we need "doubled" precision for 1/e
For 12 digits, 1/e ≈ 0.367879441171 + 0.442321595524E-12

To reduce iterations, we also need good starting estimate.
Lua e^W code (0, -1 branch): https://www.hpmuseum.org/forum/thread-19...#pid167919

That's pretty much what I feared. I imagine that the best compromise would be to use extended precision reals, but SysRPL is beyond my skill set.

The estimation code I used for branch 0 is amazingly accurate over the entire numerical range representable by Saturn-based HP's. Also, the recurrence that I used converges faster than Newton's method while not being much larger in code size. Newton's method may be somewhat more accurate near the branch point but it wouldn't be much of an improvement without double-precision numbers.

Of more interest to me is computing W(z) for branches other than 0 and -1. I have not been able to find any methods for doing so. Any insight in this area would be greatly appreciated.

[Albert Chan]

(03-21-2023 01:53 PM)John Keith Wrote:  I imagine that the best compromise would be to use extended precision reals, but SysRPL is beyond my skill set.

Just to be clear, we don't need extended precison numbers.
My Lua code only use the same 64-bits float throughout.

All we need is more precise constant 1/e, so that (y-1/e) does not suffer cancellation errors.

Quote:The estimation code I used for branch 0 is amazingly accurate over the entire numerical range representable by Saturn-based HP's. Also, the recurrence that I used converges faster than Newton's method while not being much larger in code size.

I am unfamilar with RPL, can you post actual formula for W(a) guess, and the iteration formula?

Convergence can be affected greatly by what the iteration formula look like.
Newton's method is not necessarily much slower than higher order method.

To give an example, lets do W(a = 1e10), with guess a0 = log1p(a) ≈ 23.0259
Note: nestlist is my defined function, not part of XCas

XCas> a := 1e10; a0 := log1p(a)
XCas> f := x*e^x - a
XCas> nestlist(unapply(x-f/f',x), a0, 3)                      → [23.0259, 22.1091, 21.2606, 20.568]
XCas> nestlist(unapply(x-f/(f'-f''/2*f/f'),x), a0, 3)      → [23.0259, 21.2714, 20.1788, 20.029]

This is a bad setup, with huge derivatives, f' = (x+1)*e^x, f'' = (x+2)*e^x
If a is bigger , bad setup is going to take longer to converge.

Let's try a better setup.

XCas> f := x - ln(a/x)
XCas> nestlist(unapply(x-f/f',x), a0, 3)                      → [23.0259, 20.0198, 20.0287, 20.0287]
XCas> nestlist(unapply(x-f/(f'-f''/2*f/f'),x), a0, 3)      → [23.0259, 20.0279, 20.0287, 20.0287]

Here, there is not much need for complicated halley's method. Newton is almost as good.
If a is bigger , better setup is going to converge faster! Perhaps in 1 iteration!

Quote:Of more interest to me is computing W(z) for branches other than 0 and -1. I have not been able to find any methods for doing so. Any insight in this area would be greatly appreciated.

Python mpmath lambertw(z, k=0) support all branches.

https://mpmath.org/doc/0.19/functions/po...w-function
https://github.com/mpmath/mpmath/blob/ma...ns.py#L464

[John Keith]

(03-21-2023 05:15 PM)Albert Chan Wrote:  I am unfamilar with RPL, can you post actual formula for W(a) guess, and the iteration formula?


Python mpmath lambertw(z, k=0) support all branches.

https://mpmath.org/doc/0.19/functions/po...w-function
https://github.com/mpmath/mpmath/blob/ma...ns.py#L464

Many thanks for the mpmath github reference, that solved my branch problem! Updated program will follow soon.

Here are the formulas that you requested:

Formula for global approximation of W(z) for branch 0.
From "New approximations to the principal real-valued
branch of the Lambert W-function", by Iacono and Boyd.

-1 + a*ln((1 + b*y)/(1 + c*ln(1+y))

where

y = sqrt(1+e*x)

c = (e^(1/a) - 1 - sqrt(2)/a)/(1 - e^(1/a)*ln(2))

b = sqrt(2)/a + c


With 'a' determined empirically to be 2.036, final formula is:

-1 + 2.036*ln((1 + 1.14956131*y)/(1 + 0.4549574*ln(1 + y))


Recurrence from same paper, also from "Guaranteed- and high-precision
evaluation of the Lambert W function" by Lajos Lóczi.

B'(x) = (B(x)/(1 + B(x))*(1 + ln(x/B(x))

where

B(x) is current estimate
B'(x) is new estimate

[John Keith]

(03-21-2023 06:15 AM)Gerald H Wrote:  This programme works for positive reals - I'm interested in the comparison of accuracy:

They are roughly comparable, unless x is very large. A few examples:

Code:


x = -.367
Exact:       -.932399184748
John Keith:  -.932399184762
Gerald H:    -.932399184749

x = -.36787944
Exact:       -.999920198484
John Keith:  -.99992018987
Gerald H:    -.999920176091

x = -.367879441171
Exact:       -.999998449288
John Keith:  -.99999787867
Gerald H:    -.999998020036

x = .000001
Exact:       9.99999000001E-7
John Keith:  .000000999999
Gerald H:    9.99999000001E-7

x = 999999999999
Exact:       24.4350044049
John Keith:     same
Gerald H:       same

x = 9.99999999999E499
Exact:       1144.25004187
John Keith:     same
Gerald H:    1144.25102833

Exact values were computed with Mathematica version 12.2.

As you can see, both programs lose accuracy for values very close to -1/e. Your program is very accurate over most of the real range but less accurate at MAXR. This is a known problem with Newton's method.

[Albert Chan]

(03-22-2023 07:30 PM)John Keith Wrote:  y = sqrt(1+e*x)
...
With 'a' determined empirically to be 2.036, final formula is:

-1 + 2.036*ln((1 + 1.14956131*y)/(1 + 0.4549574*ln(1 + y))

Recurrence from same paper, also from "Guaranteed- and high-precision
evaluation of the Lambert W function" by Lajos Lóczi.

B'(x) = (B(x)/(1 + B(x))*(1 + ln(x/B(x))

There is a flaw with Iacono and Boyd W guess formula.
If x is tiny enough, below 6e-10, it could return a negative guess.
This cause Newton formula to take a log of negative number, which crash it.

[Albert Chan]

Here is my implementation of accurate x = W(a), both branches.

x * e^x = a
ln(x) + x = ln(a)      → f = x + ln(x/a) = 0

x - f/f' = x - (x + ln(x/a)) * x/(1+x) = (1 - ln(x/a)) * x/(1+x)

W(-1/e) = -1 --> a ≈ -1/e, x ≈ -1 --> (1+x) is tiny

Let 1/e = r + err

1 - ln(x/a) = - ln(x*(r+err)/a) = - log1p((x*err-(a+r)+r*(1+x))/a)

x + ln(x/a) = (x+1) - (1 - ln(x/a))

Terms inside log1p sorted in increasing size, r*(1+x) the biggest.
With this, and some minor changes to my expW function, I get this:

Code:

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

lua> W(-0.367, 0, true)
-0.9323991847479226
-0.9323991847479282
lua> W(-0.367, -1, true)
-1.0707918867680617
-1.0707918867680521

lua> W(1e10, 0, true)
18.111645253927524
20.02372402639066
20.028685384073526
20.02868541330495
20.028685413304952

Update 1: added special case W₀(0) = 0, W_-1(0) = -∞ (should it be NaN?)

Update 2: removed guard "if h == err then return -1 end", see next post.

Update 3: improve W guess with a Newton's step, based from y = e^W
y = (y+a) / (log(y)+1)      → x = a/y = a/(y+a) * (log(y)+1)

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

[John Keith]

(03-23-2023 12:16 AM)Albert Chan Wrote:  There is a flaw with Iacono and Boyd W guess formula.
If x is tiny enough, below 6e-10, it could return a negative guess.
This cause Newton formula to take a log of negative number, which crash it.

That explains the problem that I saw when |z| was < 10^8. Lines 3 and 7 of my first program are a patch around the flaw. The program doesn't crash without the patch but it returns a complex result in a different branch.

There is one line in your program that concerns me:

if h == err then return -1 end

It seems that the program would return -1 if a is the closest machine representation of -1/e, which would be -.367879441171 on the HP calculators. Instead, we would want the program to return -.999998449288 in branch 0, or -1.00000155071 in branch -1, which are the correct 12-digit representations of W(-.367879441171).

My apologies if I am reading your program incorrectly.

Also, there is another line whose meaning is unclear to me:

x = sqrt(2*h/r) * (x and -1 or 1)

In RPL this statement always returns 1.

[Albert Chan]

(03-23-2023 05:53 PM)John Keith Wrote:  There is one line in your program that concerns me:

if h == err then return -1 end

It seems that the program would return -1 if a is the closest machine representation of -1/e, which would be -.367879441171 on the HP calculators. Instead, we would want the program to return -.999998449288 in branch 0, or -1.00000155071 in branch -1, which are the correct 12-digit representations of W(-.367879441171).

You are right.

Guard should be if (h == 0) ... but that would never happen.
Guard is not needed. I'll remove it.

Quote:x = sqrt(2*h/r) * (x and -1 or 1)

Last term is Lua idiom for ternary operators

It says that pick -1 if x evaluated true (neither nil nor false), 1 otherwise.
Signs are the result of solving W guess quadratic equation, 2 roots for W 2 branches

[Albert Chan]

(03-23-2023 02:56 AM)Albert Chan Wrote:  Here is my implementation of accurate x = W(a), both branches.

x * e^x = a
ln(x) + x = ln(a)      → f = x + ln(x/a) = 0

x - f/f' = x - (x + ln(x/a)) * x/(1+x) = (1 - ln(x/a)) * x/(1+x)

I just realized this is equivalent to solving for y = e^W(a)
With "same" guess, x*y=a, Newton convergence rate are identical.

x = (1 - ln(x/a)) * x/(1+x)
a/y = (1 + ln(y)) * 1/(1+y/a)

y = (y+a) / (1 + ln(y))

lua> expW(1e99, nil, true)
2.5e+098      4
5.492824331831047e+096      182.05570387623288
4.493790313227532e+096      222.52929716290643
4.493356750520384e+096      222.55076895111614
4.493356750426822e+096      222.55076895575016
4.493356750426821e+096      222.55076895575021

lua> W(1e99, nil, true)
182.05570387623249
222.52929716290646
222.55076895111617
222.55076895575016
222.5507689557502

Note: W guess used 1 Newton step. Both code converged the same way.
Because solving for y = e^W(a) is easier, W code is not recommended.

Bonus: with y, there are 2 ways to recover W(a) = a/y or log(y)
Depends on size of y (use log if y is big), we can make W(a) estimate better.

Example, above last 2 y's differ by 1 ULP, but the log's are the same.
(Not exactly. Here, y 1 ULP error translated to 0.007 ULP error in x)

lua> log(4.493356750426822e+096)
222.5507689557502
lua> log(4.493356750426821e+096)
222.5507689557502

[John Keith]

Hi Albert,

I am still having problems with your Lua program. This line, h = x/(1+x) * s(x) returns zero or a very small number (usually 10^-12 or less) and thus no correction occurs in the line x = x - h. This may be caused by the difference between the 64-bit binary floats used by Lua and the 12-digit BCD floats used by HP calculators. It is also possible that I made an error in translating the program although I have checked carefully.

I was careful to sum the values in this line s = function(x) return (x+1) + log1p((x*err-(a+r)+r*(1+x))/a) in the proper order as suggested in your post.

I have also seen several mentions of using Taylor series rather than iteration for values close to -1/e. This may also be worth exploring.

[Albert Chan]

(03-27-2023 04:45 PM)John Keith Wrote:  This line, h = x/(1+x) * s(x) returns zero or a very small number (usually 10^-12 or less)

Hi, John Keith

I assumed you meant correction around branch-point is small. That's normal.
My W guess formula is pretty accurate for a ≈ -1/e, both branches.

Solve s(x) = x + ln(x/a) = 0 (or its messy version) is probably a bad idea.
My previous post suggested it is better to solve for y = e^W(a), shown below.

lua> a = -0.3678 -- calculate W(a), 0 branch, by "hand"
lua> r, err = 0.36787944117144233, -1.2428753672788363E-17
lua> h = a + r + err
lua> x = sqrt(2*h/r) -- rough guess for e*h = (1+x)*log1p(x) - x
lua> x = x * sqrt(1+x/3) -- better guess
lua> y = r + r*x -- = e^W(a), if x is correct
lua> y, a/y, log(y)
0.3755511062373703      -0.9793607152032429      -0.9793607152084147

W(a) 2 ways matched does not imply they are correct, only that we are close.
Interestingly, true W(a) may not be a "mean" of the two ways.

lua> y = (y+a) / log1p((y-r-err)/r) -- Newton's step
lua> y, a/y, log(y)
0.37555110633147754    -0.9793607149578304      -0.9793607149578304

y's had doubled in accuracy, to full precision = e^W(a)

Quote:I have also seen several mentions of using Taylor series rather than iteration for values close to -1/e.

XCas> H := (1+x)*ln(1+x) - x
XCas> x0 := sqrt(2*H)               /* rough guess for above formula */
XCas> r0 := revert(series(x0,0,8))

\(\displaystyle x
+\frac{x^2}{6}
-\frac{x^3}{72}
+\frac{x^4}{270}
-\frac{23 x^5}{17280}
+\frac{19 x^6}{34020}
-\frac{11237 x^7}{43545600}
+\frac{13 x^8}{102060}\)

XCas> rs := series((r0/x)^2)      /* correction, inside square root */

\(\displaystyle
1+\frac{x}{3}
+\frac{x^3}{360}
-\frac{x^4}{810}
+\frac{23 x^5}{40320}
+ O\!\left(x^6\right)\)

This explained why my W guess formula is good (no x^2 terms, others tiny, with opposite sign)

Let's numerically check if series is good, with pade approximation matching to x^4
(replace 57 by 4032/71 ≈ 56.8, pade matched to x^5, but make no difference here)

lua> x = sqrt(2*h/r) -- rough guess for e*h = (1+x)*log1p(x) - x
lua> x = x * sqrt(1+x/3 + x^3/(360+160*x*(1-x/57)))
lua> y = r + r*x -- = e^W(a), if x is correct
lua> y, a/y, log(y)
0.3755511063314775      -0.9793607149578305      -0.9793607149578305

For comparison, my W code use simple guess (no bolded mess) + 1 Newton step.
Newton's step is very cheap; we had to convert x to W(a) guess anyway.

W(a) = a/y = a * (ln(y)+1)/(y+a), with y = r+r*x, we have:

Code:

x = a * log1p(x)/(r*x+h) -- Newton step

lua> W(-0.3678, nil, true) -- h=0 on the first try.
-0.9793607149578305
-0.9793607149578305

[Albert Chan]

(03-27-2023 08:57 PM)Albert Chan Wrote:  lua> h = a + r + err
lua> x = sqrt(2*h/r) -- rough guess for e*h = (1+x)*log1p(x) - x
lua> x = x * sqrt(1+x/3) -- better guess
lua> y = r + r*x -- = e^W(a), if x is correct

Technically, it should be y = (1+x)/e = (1+x)*(r+err)
Numerically, err < 1/2 ULP of r, thus y = r + r*x

e*h = (e*y) * log(e*y) - (e*y-1)
h = y * (log(y)+1) - y + 1/e
h - 1/e = y*log(y) = a

---

A novel way is to directly solve for x, using log1p_sub(x), then convert to W(a) or e^W(a)
Except for getting accurate h, err = 1/e - float(1/e) is not used anymore.

Let H=e*h, L = log(1+x) - x, accurately computed

f = (1+x)*log(1+x) - x - H
f' = log(1+x) + (1+x)/(1+x) - 1 = log(1+x)

x - f/f' = (H-L) / (x+L)

Note that numerator is really an addition, with both H, -L positive.

lua> x -- better guess, from above quote
0.020853747742736108
lua> H = h/r
lua> L = log1p_sub(x)
lua> x = (H-L) / (x+L) -- newton step
lua> x -- converged to full precision
0.020853747998545925
lua> r+r*x, log1p(x)-1 -- = e^W(a), W(a)
0.3755511063314775      -0.9793607149578305

[Albert Chan]

(03-27-2023 08:57 PM)Albert Chan Wrote:  lua> x = sqrt(2*h/r) -- rough guess for e*h = (1+x)*log1p(x) - x
lua> x = x * sqrt(1+x/3) -- better guess
lua> y = r + r*x -- = e^W(a), if x is correct
lua> y, a/y, log(y)
0.3755511062373703      -0.9793607152032429      -0.9793607152084147

W(a) 2 ways matched does not imply they are correct, only that we are close.
Interestingly, true W(a) may not be a "mean" of the two ways.

I was curious why the 2 W's do not bracket true W(a)
Doing a bit of error analysis, everything becomes clear.

Ignoring error of division, we have: relerr(x = a/y) ≈ - relerr(y)

ln(y*(1+ε)) ≈ ln(y) + ε
ln(y*(1+ε)) / ln(y) ≈ 1 + ε/ln(y)     relerr(ln(y)) ≈ relerr(y) / ln(y)

2 W's don't bracket true W(a) is because ln(y) ≈ x is negative.
2 W's have errors of the same sign.

relerr(x) + x*relerr(ln(y)) ≈ -relerr(y) + relerr(y) ≈ 0

lua> x, lny = -0.9793607152032429, -0.9793607152084147
lua> (x + x*lny) / (1+x)
-0.9793607149578298

Simply extrapolate for 0 error, we get good W(-0.3678) estimate.
Actually, extrapolate for 0 error is equivalent to a Newton's step!

[John Keith]

(03-31-2023 04:06 PM)Albert Chan Wrote:  lua> x, lny = -0.9793607152032429, -0.9793607152084147
lua> (x + x*lny) / (1+x)
-0.9793607149578298

Simply extrapolate for 0 error, we get good W(-0.3678) estimate.
Actually, extrapolate for 0 error is equivalent to a Newton's step!

Thanks for your help and insights, this is all quite interesting if a bit beyond my level of knowledge.

Unfortunately, I am still not able to achieve that level of accuracy with the 12-digit precision of HP calculators.

For an input of -.3678, your program for e^W(x) returns .375551106194 and -.979360715069

Applying your correction above returns -.979360714941 which is a noticeable improvement but still off by 16 ULP's. My program from post #1 returns -.979360714903 which is a bit worse (off by 54 ULP's) but with no correction applied.

What I am hoping for is a way to achieve 11-digit accuracy for values very close to -1/e, such as -.36787944117. I still don't think this will be possible without extended precision but I am hoping to be pleasantly surprised.

I am close to completing a new version of my program which covers the entire complex plane in all branches. It seems to be accurate to within 2 ULP's except for a circle of radius ~.001 around -1/e. It is only that area that i am concerned about.

[Albert Chan]

(03-31-2023 06:15 PM)John Keith Wrote:  Unfortunately, I am still not able to achieve that level of accuracy with the 12-digit precision of HP calculators.

For an input of -.3678, your program for e^W(x) returns .375551106194 and -.979360715069

My Lua expW code translated to HP71B. It seems to work fine.

Code:

10 INPUT "a, k? ";A,K @ K=K+K+1
20 IF A<=-.1 THEN
30 R=.367879441171 @ R2=4.42321595524E-13
40 Y=SQRT(2*R*(A+R+R2))*K @ Y=R+Y*SQRT(1+Y/(3*R))
50 REPEAT @ H=Y-(Y+A)/LOGP1((Y-R-R2)/R) @ Y=Y-H @ UNTIL Y=Y+H*.0001
60 ELSE
70 Y=(K=1)+K*A/4
80 REPEAT @ H=Y-(Y+A)/(LOG(Y)+1) @ Y=Y-H @ UNTIL Y=Y+H*.0001
90 END IF
100 PRINT "e^W, W =";Y;A/Y

>run
a, k? -.3678, 0
e^W, W = .37555110633 -.979360714962
>run
a, k? -.3678, -1
e^W, W = .360260737301 -1.02092723941

>run
a, k? -1e-6, 0
e^W, W = .999998999999 -.000001000001
>run
a, k? -1e-6, -1
e^W, W = 6.01449171278E-8 -16.6265089014

[Albert Chan]

(03-27-2023 11:30 PM)Albert Chan Wrote:  Let H=e*h, L = log(1+x) - x, accurately computed

f = (1+x)*log(1+x) - x - H
f' = log(1+x) + (1+x)/(1+x) - 1 = log(1+x)

x - f/f' = (H-L) / (x+L)

Note that numerator is really an addition, with both H, -L positive.
...
lua> r+r*x, log1p(x)-1 -- = e^W(a), W(a)

When we are close to branch point, x is tiny. Even rough x estimate suffice.
When we are far away, log1p(x)-x does not suffer catastrophic cancellation.

We may not need accurate log1p(x)-x after all.

Here, I added code to roughly solve above f=0 for x, then convert to W's
Old code with iterations of y remained (2nd e^W, W), for comparison.

>40 H=(A+R+R2)/R @ X=SQRT(2*H)*K @ X=X*SQRT(1+X/3) @ Y=R+R*X
>42 REPEAT @ Z=LOGP1(X) @ Z=X-(X-Z+H)/Z @ X=X-Z @ UNTIL X=X+Z*.000001
>44 PRINT "e^W, W =";R+R*X;LOGP1(X)-1

>run
a, k? -.367879441171,0
e^W, W = .367880011645 -.999998449291 ! true W = -.999998449288
e^W, W = .367880011646 -.999998449291
>run
a, k? -.36787944117,0
e^W, W = .367880471317 -.999997199776 ! true W = -.999997199775
e^W, W = .367880471317 -.999997199778
>run
a, k? -.3678794411,0
e^W, W = .367886691321 -.999980292244 ! true W = -.999980292245
e^W, W = .36788669132   -.999980292247
>run
a, k? -.367879441,0
e^W, W = .367890672444 -.999969470701 ! true W = -.999969470701
e^W, W = .367890672444 -.999969470702

...

>run
a, k? -.3678,0
e^W, W = .375551106332 -.979360714956 ! true W = -.979360714958
e^W, W = .37555110633   -.979360714962

[Albert Chan]

(03-31-2023 10:07 PM)Albert Chan Wrote:  We may not need accurate log1p(x)-x after all.

With accurate log1p(x)-x, we can get W around branch point to within 1 ULP or less.
With below FNL(X), previous post examples all match W exactly, no ULP errors.

200 DEF FNL(X) ! = ln(1+X) - X, but more accurate
210 IF ABS(X)>=.4 THEN X=X/(SQRT(1+X)+1) @ FNL=FNL(X)*2-X*X @ GOTO 250
220 X2=X/(X+2) @ X4=X2*X2
230 X4=X4*(5005-X4*(5082-X4*969))/(15015-X4*(24255-X4*(11025-X4*1225)))
240 FNL=X2*(X4+X4-X)
250 END DEF

line 230: X4 = pade approx. of atanh(X2)/X2 - 1 ≈ ln(1+X)/(2*X2) - 1
line 240: FNL = X2*(X4+X4-X) ≈ ln(1+X) - 2*X2 - X*X2 = ln(1+X) - X

42 REPEAT @ L=FNL(X) @ Z=X-(H-L)/(X+L) @ X=X-Z @ UNTIL X=X+Z*.0001
44 PRINT "e^W, W =";R+R*X;LOGP1(X)-1

Away from branch point is a different story. With X getting big, FNL(X) won't help much.