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
 Albert Chan Senior Member Posts: 2,699 Joined: Jul 2018
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

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
10-03-2020, 05:29 PM
Post: #42
 Albert Chan Senior Member Posts: 2,699 Joined: Jul 2018
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
10-03-2020, 07:56 PM
Post: #43
 lyuka Junior Member Posts: 48 Joined: May 2014
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.
10-05-2020, 08:03 AM
Post: #44
 Werner Senior Member Posts: 891 Joined: Dec 2013
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
10-05-2020, 03:20 PM
Post: #45
 Albert Chan Senior Member Posts: 2,699 Joined: Jul 2018
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
10-05-2020, 06:09 PM (This post was last modified: 10-06-2020 09:52 AM by lyuka.)
Post: #46
 lyuka Junior Member Posts: 48 Joined: May 2014
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.
10-06-2020, 06:16 AM (This post was last modified: 10-06-2020 06:16 AM by Werner.)
Post: #47
 Werner Senior Member Posts: 891 Joined: Dec 2013
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
10-21-2020, 03:05 PM (This post was last modified: 03-30-2021 10:48 AM by Albert Chan.)
Post: #48
 Albert Chan Senior Member Posts: 2,699 Joined: Jul 2018
RE: Lambert W Function (hp-42s)
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^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"
11-09-2020, 08:30 AM
Post: #49
 lyuka Junior Member Posts: 48 Joined: May 2014
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 :-)
08-24-2022, 09:48 PM (This post was last modified: 08-25-2022 01:57 AM by Bill Triplett.)
Post: #50
 Bill Triplett Junior Member Posts: 24 Joined: Dec 2013
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
12-29-2022, 04:17 AM (This post was last modified: 12-29-2022 09:08 PM by Albert Chan.)
Post: #51
 Albert Chan Senior Member Posts: 2,699 Joined: Jul 2018
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
12-29-2022, 12:21 PM
Post: #52
 Albert Chan Senior Member Posts: 2,699 Joined: Jul 2018
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
12-29-2022, 08:01 PM (This post was last modified: 12-30-2022 01:24 AM by Albert Chan.)
Post: #53
 Albert Chan Senior Member Posts: 2,699 Joined: Jul 2018
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
 « Next Oldest | Next Newest »

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