Post Reply 
Lambert W Function (hp-42s)
10-03-2020, 04:31 PM (This post was last modified: 10-03-2020 04:51 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
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
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
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
Post Reply 




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