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 Code: 00 { 56-Byte Prgm } 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
|
|||
|
|||
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
|
|||
|
|||
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
|
|||
|
|||
RE: Lambert W Function (hp-42s)
This time, I'll stick to what I'm good at:
Code: 00 { 53-Byte Prgm } Cheers, Werner 41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE |
|||
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" 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" |
|||
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:
[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
|
|||
|
|||
RE: Lambert W Function (hp-42s)
I can make it shorter still, but you won't like it:
Code: 00 { 51-Byte Prgm } 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
|
|||
|
|||
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
|
|||
|
|||
RE: Lambert W Function (hp-42s) | |||
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 } |
|||
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 } 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
|
|||
|
|||
RE: Lambert W Function (hp-42s) | |||
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 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 } |
|||
« Next Oldest | Next Newest »
|
User(s) browsing this thread: 2 Guest(s)