Lambert W Function (hp-42s) - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html) +--- Forum: General Forum (/forum-4.html) +--- Thread: Lambert W Function (hp-42s) (/thread-15004.html) |
RE: Lambert W Function (hp-42s) - Albert Chan - 10-03-2020 04:31 PM 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 RE: Lambert W Function (hp-42s) - Albert Chan - 10-03-2020 05:29 PM (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 RE: Lambert W Function (hp-42s) - lyuka - 10-03-2020 07:56 PM 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. RE: Lambert W Function (hp-42s) - Werner - 10-05-2020 08:03 AM This time, I'll stick to what I'm good at: Code: 00 { 53-Byte Prgm } Cheers, Werner RE: Lambert W Function (hp-42s) - Albert Chan - 10-05-2020 03:20 PM 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" RE: Lambert W Function (hp-42s) - lyuka - 10-05-2020 06:09 PM (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. RE: Lambert W Function (hp-42s) - Werner - 10-06-2020 06:16 AM I can make it shorter still, but you won't like it: Code: 00 { 51-Byte Prgm } Cheers, Werner RE: Lambert W Function (hp-42s) - Albert Chan - 10-21-2020 03:05 PM 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" RE: Lambert W Function (hp-42s) - lyuka - 11-09-2020 08:30 AM 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 :-) RE: Lambert W Function (hp-42s) - Bill Triplett - 08-24-2022 09:48 PM 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 } RE: Lambert W Function (hp-42s) - Albert Chan - 12-29-2022 04:17 AM 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 RE: Lambert W Function (hp-42s) - Albert Chan - 12-29-2022 12:21 PM (12-29-2022 04:17 AM)Albert Chan Wrote: If y' is big, this does not matter, our test implied relative error test too. Later may mean infinte loops ... so probably not OK RE: Lambert W Function (hp-42s) - Albert Chan - 12-29-2022 08:01 PM (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 } |