Post Reply 
Lambert W Function (hp-42s)
10-03-2020, 04:31 PM (This post was last modified: 12-30-2022 10:41 PM by Albert Chan.)
Post: #41
RE: Lambert W Function (hp-42s)
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 Smile

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
Find all posts by this user
Quote this message in a reply
10-03-2020, 05:29 PM
Post: #42
RE: Lambert W Function (hp-42s)
(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
Find all posts by this user
Quote this message in a reply
10-03-2020, 07:56 PM
Post: #43
RE: Lambert W Function (hp-42s)
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.
Find all posts by this user
Quote this message in a reply
10-05-2020, 08:03 AM
Post: #44
RE: Lambert W Function (hp-42s)
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

41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE
Find all posts by this user
Quote this message in a reply
10-05-2020, 03:20 PM
Post: #45
RE: Lambert W Function (hp-42s)
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
Find all posts by this user
Quote this message in a reply
10-05-2020, 06:09 PM (This post was last modified: 10-06-2020 09:52 AM by lyuka.)
Post: #46
RE: Lambert W Function (hp-42s)
(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.
Find all posts by this user
Quote this message in a reply
10-06-2020, 06:16 AM (This post was last modified: 10-06-2020 06:16 AM by Werner.)
Post: #47
RE: Lambert W Function (hp-42s)
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

41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE
Find all posts by this user
Quote this message in a reply
10-21-2020, 03:05 PM (This post was last modified: 03-30-2021 10:48 AM by Albert Chan.)
Post: #48
RE: Lambert W Function (hp-42s)
Hi, lyuka

I discovered a feature of your code. Smile
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^W0(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"
Find all posts by this user
Quote this message in a reply
11-09-2020, 08:30 AM
Post: #49
RE: Lambert W Function (hp-42s)
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 :-)
Find all posts by this user
Quote this message in a reply
08-24-2022, 09:48 PM (This post was last modified: 08-25-2022 01:57 AM by Bill Triplett.)
Post: #50
RE: Lambert W Function (hp-42s)
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
Find all posts by this user
Quote this message in a reply
12-29-2022, 04:17 AM (This post was last modified: 12-29-2022 09:08 PM by Albert Chan.)
Post: #51
RE: Lambert W Function (hp-42s)
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"    → W0(-ε) = 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
Find all posts by this user
Quote this message in a reply
12-29-2022, 12:21 PM
Post: #52
RE: Lambert W Function (hp-42s)
(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 Big Grin
Find all posts by this user
Quote this message in a reply
12-29-2022, 08:01 PM (This post was last modified: 12-30-2022 01:24 AM by Albert Chan.)
Post: #53
RE: Lambert W Function (hp-42s)
(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
Find all posts by this user
Quote this message in a reply
Post Reply 




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