Lambert W Function (hp-42s)

[Albert Chan]

We were testing special case x = -r = -1/e
I was curious, is it the same as guess(x) = -x ?

guess(x) = r + 0.3*(x+r) + √(2r*(x+r)) = -x   → √(2r*(x+r)) = -1.3*(x+r)

Let z = x+r, we have √(2/e*z) = -1.3*z           → z = 0 → x = -r

With this, we can simplify the code greatly, moving the test after guess(x).
And, turning repeat until loop into while do loop, remove the test altogether !

Bonus: now it checked both special case, x = -1/e, -1/e + 0j

Code:

00 { 56-Byte Prgm }
01▶LBL "eW"
02 -1
03 E↑X          # r     x
04 ENTER
05 RCL+ ST Z    # x+r   r       x
06 ENTER
07 STO+ ST X    # 2x+2r x+r     r   x
08 RCL× ST Z
09 SQRT         # sqrt  x+r     r   x
10 STO+ ST Z
11 R↓           # x+r   sqrt+r  x   sqrt
12 0.3
13 ×
14 +            # y     x       x   x
15 X<>Y
16 +/-
17 X<>Y         # y     -x      x   x
18▶LBL 01
19 X=Y?         # Convergence check
20 RTN
21 STO ST Y
22 LN
23 1
24 +
25 R↑           
26 RCL+ ST Z
27 X<>Y         
28 ÷            # Newton-Raphson method
29 -            
30 STO× ST X    
31 LASTX        
32 STO+ ST Y    # y'    y'+dy^2 x   x
33 GTO 01
34 END

UPDATE:
There is a flaw in termination test.
I keep this for historical record. Please use corrected version

[Albert Chan]

(10-03-2020 04:31 PM)Albert Chan Wrote:  Let z = x+r, we have √(2/e*z) = -1.3*z           → z = 0 → x = -r

I tried to solve by hand. Can someone check if this is correct ?

√(2/e*z) = -1.3*z
√z * (√(2/e) + 1.3*√z) = 0

√z = 0, or -√(2/e)/1.3 ≈ -0.6598

Because of square root branch cut, Re(√z) ≥ 0, thus we are left with z = 0

[lyuka]

Hi Albert,

Thanks to you, the eW program is now so sophisticated.

I have confirmed that, for complex number x, the guess y become -x only when x = 0.

[Werner]

This time, I'll stick to what I'm good at:

Code:

00 { 53-Byte Prgm }
01▸LBL "eW"
02 0.3 @        L       X       Y       Z       T
03 -1
04 E^X @                r       .3      x
05 RCL+ ST Z @  r       r+x     .3      x
06 STO× ST Y @                  .3(r+x)
07 STO+ ST X @          2r+2x
08 LASTX @              r       2r+2x   .3(r+x) x
09 STO+ ST Z @          r       2r+2x   r+"     x
10 ×
11 SQRT
12 +
13 X<>Y
14 +/-
15 X<>Y @               y       -x      x       x
16▸LBL 01
17 X=Y? @ Convergence check
18 RTN
19 STO ST Y
20 LN
21 1
22 +
23 R^
24 RCL+ ST Z
25 X<>Y
26 ÷ @ Newton-Raphson method
27 -
28 STO× ST X
29 LASTX
30 STO+ ST Y @ y'    y'+dy^2 x   x
31 GTO 01
32 END

Cheers, Werner

[Albert Chan]

Hi Werner

Nice idea putting constant 0.3 up front !

Tracking LASTX is hard, it might be a good idea to use it up ASAP, like this:
(I am leaving out the LBL01 ... GTO 01 loops)

Code:

01▸LBL "eW"
02 0.3 @        L       X       Y       Z           T
03 -1
04 E^X @                r       .3      x
05 RCL+ ST Z @  r       r+x     .3      x
06 STO× ST Y @                  .3(r+x)
08 LASTX @              r       r+x     .3(r+x)     x
09 STO+ ST Z @         
09 STO+ ST X @          2r      r+x     .3(r+x)+r   x
10 ×
11 SQRT
12 +
13 X<>Y
14 +/-
15 X<>Y @               y       -x      x       x

Here is a version that does not use LASTX to build guess, (same steps and bytes count)
Rewrite guess(x) = √(2r(x+r)) + 1.3(x+r) - x

This is not as accurate when |x| is huge, but for guess, it does not matter.
The up side is that code is easier to follow.

Code:

01▸LBL "eW"
02 1.3
03 -1
04 E^X
05 ENTER
07 RCL+ ST T @      #   r+x     r      1.3      x
08 STO× ST Z @
09 STO+ ST X @      #   2(r+x)  r   1.3(r+x)    x
09 ×
10 SQRT       
11 +                #   y+x     x       x       x
12 X<>Y
13 +/-              #   -x      y+x     x       x
14 STO ST Z
15 +                #    y      -x      x       x

[lyuka]

(10-05-2020 03:20 PM)Albert Chan Wrote:  This is not as accurate when |x| is huge, but for guess, it does not matter.

BTW, In my function of the eW(x) for real x written in C, the following code is used for the guess.

Code:


#include <math.h>
r = 1.0 / M_E;
y = x < 7.9488901779546342
        ? r + sqrt((r + r) * (x + r)) + 0.3 * (x + r)
        : 1.1 * x / log(0.5 * x) - 1.0; /* Initial value */

[attachment=8779]
This guess y is accurate within +/- 0.09 for -1/e < x < 1E128, and it may reduce one or two itteration for the Newton-Raphson method, however, it may not so useful as it requires extra condition branch and one log() calculation.

[Werner]

I can make it shorter still, but you won't like it:

Code:

00 { 51-Byte Prgm }
01▸LBL "eW"
02 +/-
03 1.3 @        L       X       Y       Z       T
04 -1
05 E↑X
06 ENTER @              r       r       1.3     -x
07 RCL- ST T @          r+x
08 STO× ST Z @                          1.3(r+x)
09 STO+ ST X @          2r+2x
10 ×
11 SQRT
12 +
13 + @               y       -x      -x       -x
14▸LBL 01
15 X=Y? @ Convergence check
16 RTN
17 STO ST Y
18 LN
19 1
20 +
21 RCL ST Y
22 RCL- ST T
23 X<>Y
24 ÷ @ Newton-Raphson method
25 -
26 STO× ST X
27 LASTX
28 STO+ ST Y @ y'    y'+dy^2 -x   -x
29 GTO 01
30 END

Cheers, Werner

[Albert Chan]

Hi, lyuka

I discovered a feature of your code.
Because we turned repeat-until to while-do loop, it is restartable with another guess.

Example, to calculate eW(-log(2)/2), both branch 0 and -1:

2 LN 2 +/- /   → a = -0.3465735902799726547086160607290883
XEQ "eW"       → e^W₀(a) = 0.5
− 0.3 R/S      → e^W_-1(a) = 0.25

For e^W_-1(a), -1/e < a < 0, any positive guess below -a work.

Update: a good guess for e^W_-1(a), -1/e < a < 0, is 2a²
For above case, 2a² already gives 0.24022

Just replace "− new-guess R/S", with this: "− − + × R/S"

[lyuka]

Hi, Albert
I'm sorry I overlooked it again.
I was unaware of this feature.

(10-21-2020 03:05 PM)Albert Chan Wrote:  Just replace "− new-guess R/S", with this: "− − + × R/S"

This is really cool :-)

[Bill Triplett]

This work is beautiful.

In case a user does not care so much about minimum program space, here is a version based on Albert Chan's post from 10-06-2020, 01:20 AM. This version outputs W(x) instead of e^W(x), and empties the stack:

Code:

00 { 64-Byte Prgm }
01▸LBL "cW"
02 0.3
03 -1
04 E^X
05 RCL+ ST Z
06 STO× ST Y
07 STO+ ST X
08 LASTX
09 STO+ ST Z
10 ×
11 SQRT
12 +
13 X<>Y
14 +/-
15 X<>Y
16▸LBL 01
17 X=Y?
18 GTO 02
19 STO ST Y
20 LN
21 1
22 +
23 R^
24 RCL+ ST Z
25 X<>Y
26 ÷
27 -
28 STO× ST X
29 LASTX
30 STO+ ST Y
31 GTO 01
32 LBL 02
33 LN
34 X<>Y
35 CLX
36 STO ST Z
37 STO ST T
38 X<>Y
39 RTN
40 END

[Albert Chan]

I discovered a bug, for termination condition.

We used y' == y' + dy^2 termination test, but it is wrong
Assuming Newton's method have quadratic convergence, it should be:

1 == 1 + (relative error)^2
1 == 1 + (dy/y')^2
y'^2 == y'^2 + dy^2

If y' is big, this does not matter, our test implied relative error test too.
It may terminate later than needed, but that's OK.

Problem is if we use this for e^W_-1(-ε), which results in y' even smaller than ε
Now, the flawed termination test quit early.

Example, for ε = 1E-99:

-1e-99 XEQ "eW"    → W₀(-ε) = 1.
- - + * R/S        → 2.201582623716450630677604534034768e-102

This quit too early (not even 1 correct digit!), we use this guess for next iteration:

CLX + R/S          → 4.281027759337267571527542476904921e-102

We now have 3 correct digits. It will take a few more iterations for convergence.

We use this formula, which is very accurate close to x = -1/e

e^W(x) ≈ 1/e + (x+1/e)/3 + sqrt ((2/e)*(x+1/e))

I am too lazy to code relative error test, and just hard coded right eps for the job.
If y and y' matched 17 digits, this implied y' converged to 34 digits.
Worse case, we have 17+ correct digits.

Code:

00 { 59-Byte Prgm }
01▸LBL "eW"
02 3
03 1/X
04 -1
05 E↑X
06 RCL+ ST Z
07 STO× ST Y
08 STO+ ST X
09 LASTX
10 STO+ ST Z
11 ×
12 SQRT
13 +
14 X<>Y
15 +/-
16 X<>Y     @   y    -x    x    x

17▸LBL 01
18 X=Y?
19 RTN
20 STO ST Y @   y     y     x    x
21 LN
22 1
23 +
24 R↑       @   x  LN(y)+1  y    x
25 RCL+ ST Z
26 X<>Y
27 ÷        @   y'    y     x    x
28 -
29 LASTX
30 X<>Y     @   dy   y'     x    x
31 1ᴇ-17
32 ×
33 X<>Y     @   y'   eps    x    x
34 +
35 LASTX    @   y'  y'+eps  x    x
36 GTO 01
37 END

With this updated eW, we get (all digits correct)

e^W_-1(-1E-99) = 4.284330166764374654650589938202028e-102
→ W_-1(-1E-99) = -233.4087152660372939791972288026527

[Albert Chan]

(12-29-2022 04:17 AM)Albert Chan Wrote:  If y' is big, this does not matter, our test implied relative error test too.
It may terminate later than needed, but that's OK.

Later may mean infinte loops ... so probably not OK

[Albert Chan]

(10-05-2020 03:20 PM)Albert Chan Wrote:  Tracking LASTX is hard, it might be a good idea to use it up ASAP, like this ...

(12-29-2022 04:17 AM)Albert Chan Wrote:  We use this formula, which is very accurate close to x = -1/e

e^W(x) ≈ 1/e + (x+1/e)/3 + sqrt ((2/e)*(x+1/e))

I am too lazy to code relative error test, and just hard coded right eps for the job.
If y and y' matched 17 digits, this implied y' converged to 34 digits.
Worse case, we have 17+ correct digits.

This updated version fixed both issues at the same time.
Termination test is now: y' == y' * (1 + relative_error^2), independent of hardware.

Code:

00 { 56-Byte Prgm }
01▸LBL "eW"
02 3            @   X        Y        Z        T
03 1/X
04 -1
05 E^X @        @   r       1/3       x
06 RCL+ ST Z    @  r+x
07 STO× ST Y    @       (r+x)/3
08 LASTX        @   r       r+x    (r+x)/3     x
09 STO+ ST Z
10 STO+ ST X    @  2r       r+x    (r+x)/3+r   x
11 ×
12 SQRT
13 +
14 X<>Y
15 +/-
16 X<>Y         @   y       -x      x       x
17▸LBL 01
18 X=Y?
19 RTN
20 STO ST Y     @   y       y       x       x
21 LN
22 1
23 +
24 R↑
25 RCL+ ST Z
26 X<>Y         @ ln(y)+1   y+x     y       x
27 ÷
28 -
29 LASTX        @  y'       dy      x       x
30 ÷
31 STO× ST X
32 LASTX        @  y'      err^2    x       x
33 STO× ST Y
34 STO+ ST Y    @  y'  y*(1+err^2)  x       x
35 GTO 01
36 END