Lambert W Function (hp42s)

09282020, 04:06 PM
Post: #21




RE: Lambert W Function (hp42s)
Oh, I was overlooked this thread.
The Lambert W function came out when I was writing a description of the Widlar current souce a few months ago, but it is not found in ordinary scientific calculators or spreadsheet software. So I wrote an e^W calculation program for 42S, and additionally in C. http://www.finetune.co.jp/~lyuka/technot...w42s.html It's nice to be able to handle complex numbers easily with 42S, but if you try to find e^W with NewtonRaphson method, it will fail very close to 1/e. So, it is desirable that the approximation error of the initial value is asymptotic to 0 at 1/e. For that reason, I chose the following formula as an approximate expression that gives the initial value. y0 = 1 / e + sqrt ((2 / e) * (x + 1 / e)) + 0.3 * (x + 1 / e); Thanks to sqrt in this equation, in 42S we can automatically get the complex y0 from x less than 1 / e. Instead of the coefficient of 0.3 in the formula above, it can be "e  sqrt2  1" which makes zero error at x = 0, or "1 / 3" from Puiseux series, but 0.3 is quite good enough for this purpose. 

09282020, 10:01 PM
Post: #22




RE: Lambert W Function (hp42s)
(09282020 04:06 PM)lyuka Wrote: It's nice to be able to handle complex numbers easily with 42S, I use e^W(x) guess of LN(1+x), it seems OK with close to 1/e Werner improved on it by doing everything in the stack Bonus: resulting stack X = Y = eW, Z = T = x. To recover W, we can do either "LN", or "R↓ ÷" With Free42, tried x = 1/e: 1 [+/] [e^X] [+/] XEQ "eW" → 0.367879441171 It worked, but result "only" have 17 good digits. Trying it in Python, eW convergence around 1/e is bad, even with good guess. >>> from cmath import * >>> g = lambda y,a: (y+a)/(log(y)+1) >>> x = 1/e >>> y = log(1+x) # eW guess >>> y (0.45867514538708193+0j) >>> for i in range(1,101): y = g(y,x); print i, y ... 1 (0.0183829712865+0.261809735996j) 2 (0.19954449329+0.194333174479j) 3 (0.292276297365+0.112700062519j) 4 (0.332518371694+0.0601902170855j) 5 (0.350846602748+0.0310495664127j) 6 (0.359529649069+0.0157626859032j) 7 (0.363746790135+0.00794072878334j) 8 (0.365823750457+0.00398519956333j) 9 (0.36685426371+0.00199630716024j) ... 100 (0.367879435928+2.54426772287e11j) 

09292020, 09:21 AM
(This post was last modified: 09292020 09:30 AM by lyuka.)
Post: #23




RE: Lambert W Function (hp42s)
I'm sorry the word fail was wrong. It meant go bad from a convergence point of view, even if it might or might not fail.
In my opinion, it is not appropriate to use ln(x + 1) as the initial guess of e^W(x) if x is very close to 1 / e because it would require so many itterations very close to 1 / e (e.g. 1 / e + 1E11). Using Free42 with the example value of x = 1 / e + 1E11 ~= 0.3678794411614423215955237701614609, it requires 22 itterations if you start with the guess of y0 = ln(x + 1) ~= 0.4586751453712621239530755133385323. On the other hand, starting with the guess of y0 = 1 / e + sqrt ((2 / e) * (x + 1 / e)) + 0.3 * (x + 1 / e) ~= 0.3678821536620134320783140520913751, only two itterations are needed. In this case the guess is already correct upto 12 digits. 

09292020, 07:07 PM
(This post was last modified: 09292020 08:17 PM by Albert Chan.)
Post: #24




RE: Lambert W Function (hp42s)
(09282020 04:06 PM)lyuka Wrote: y0 = 1 / e + sqrt ((2 / e) * (x + 1 / e)) + 0.3 * (x + 1 / e); Very nice approximation ! However, the formula might be too good. At x = 1/e, above return *exactly* y0 = 1/e. With f(y) = y*ln(y)  x, f'(y) = ln(y) + 1, we get f'(y0) = 1 + 1 = 0 Unless we test for zero slope, Newton's formula, y = f(y)/f'(y), will crash with dividebyzero. You might want to adjust the formula a tiny bit ... 

09292020, 11:17 PM
Post: #25




RE: Lambert W Function (hp42s)
That case of x = 1 / e has already been handled by my code shown in the link in the previous post.


09302020, 07:08 AM
(This post was last modified: 09302020 07:09 AM by Werner.)
Post: #26




RE: Lambert W Function (hp42s)
(09282020 04:06 PM)lyuka Wrote: http://www.finetune.co.jp/~lyuka/technot...w42s.htmlIn that program: Code: 12 ENTER can be replaced by Code: 12 ENTER and Code: 16 X<>Y by Code: 16 STO+ ST Z saving two bytes! ;) Cheers, Werner 

09302020, 09:12 AM
Post: #27




RE: Lambert W Function (hp42s)
(09282020 10:01 PM)Albert Chan Wrote: With Free42, tried x = 1/e:That is a direct consequence of using y'=y'+(y'y)^2 as a stopping criterion.. You should do 1 extra NewtonRhapson step Werner 

09302020, 11:04 AM
Post: #28




RE: Lambert W Function (hp42s)
Thanks to Werner, I confirmed that the code was reduced by 2 bytes and updated the web page :)


09302020, 05:28 PM
Post: #29




RE: Lambert W Function (hp42s)
(09302020 09:12 AM)Werner Wrote: That is a direct consequence of using y'=y'+(y'y)^2 as a stopping criterion.. If y' = y' + (y'y)^2, it implied y has at least halfprecision. Assuming Newton's step doubled its precision, y' has about full precision. However, this does not apply on singularity. For y = e^W(x), the newton formula is: y ← y  (y*ln(y)  x) / (ln(y) + 1) Or, equivalent version (simpler, but slightly less accurate): y ← (y + x) / (ln(y) + 1) As y approach 1/e, slope (denominator) goes to zero. With limited precision, as soon as y reached half precision, ln(y) + 1 will lose half precison too. First Newton formula suffered the same issue, with questionable correction term. Both form we are hit with 0/0 issue. However, limit do exist: XCas> x0 := 1/e XCas> limit([y  (y*ln(y)x0)/(ln(y)+1), (y+x0)/(ln(y)+1)], y=1/e) → [1/e, 1/e] Let's try e^W(10^34  1/e), using lyuka "eW": .3678794411714423215955237701614608 x = 10^34  1/e = 10^34  r .3678794411714423301731626197685289 g = r + sqrt(2*r*(x+r)) + 0.3*(x+r) .3678794411714423295023829145484696 y = newton(s) w/ guess g .3678794411714423286399438752146243055... Mathematica for e^W(x) Newton's step(s) gives back only 17 good digits (as expected) Had we apply 1 more newton step (with last y) .9999999999999999785069284676275602 ln(y) .0000000000000000214930715323724398 ln(y) + 1 .3678794411714423215955237701614608 y ln(y) We have only halfprecision slope. Also, y ln(y)  x = 0. Newton steps does nothing ... 

09302020, 05:40 PM
(This post was last modified: 09302020 05:42 PM by Albert Chan.)
Post: #30




RE: Lambert W Function (hp42s)
Hi, lyuka
I extended your code to handle complex number input. Code: 00 { 57Byte Prgm } Example: e^W(1+2j) 1 Enter 2 [complex] XEQ "eW" 1.963022024795710615403657609352963 + 1.157904818620494603558662182506434i Stack X = Y = eW, Z = T = x. Confirm: x  W e^W should be tiny [LN] [×] [−] → 5.e34  1.e33i 

09302020, 07:16 PM
Post: #31




RE: Lambert W Function (hp42s)
Hi, Albert Chan.
Thanks for your complex number input extension of the code, which is elegant. I confirmed the code and updated the web page :) 

10012020, 09:37 AM
Post: #32




RE: Lambert W Function (hp42s)
nitpicking: it won't recognize (1/e,0).
We can add a code slice at the beginning that will turn a complex number with zero imaginary part into a real, as follows: Code: @ Re > Re or Code: REAL? Cheers, Werner 

10012020, 12:33 PM
(This post was last modified: 10012020 01:06 PM by Albert Chan.)
Post: #33




RE: Lambert W Function (hp42s)
(10012020 09:37 AM)Werner Wrote: nitpicking: it won't recognize (1/e,0). I don't like testing special cases, if I can avoid it. Better approach is to not returning guess y0 = 1/e, even if x ≈ 1/e One way is to take advantage of calculator internal extra precision of PI (lets call it π) Free42: [PI] [SIN] → 1.158028306006248941790250554076922e34 SIN(PI) = SIN(π  PI) ≈ π  PI Instead of adding an ε to 1/e, we get more bang for the buck if ε added to (x+1/e) Thus, another approach is to remove testing x = 1/e altogether, replacing Line 1 to 11 to this: Code: 01▸LBL "eW" To have guess returning r = 1/e, we need: y0 = r + √(2r*(x+r+ε)) + 0.3 (x+r+ε) = r → √(x+r+ε) [ √(2r) + 0.3 √(x+r+ε) ] = 0 → √(x+r+ε) = 0 or (√(2r)/0.3 ≈ 2.859213) When (x+r) approaching ε, (x+r) will suffer massive cancellations, thus x+r+ε ≠ 0 √(x+r+ε) > 0 or gone complex. 

10012020, 01:39 PM
(This post was last modified: 10012020 01:47 PM by Werner.)
Post: #34




RE: Lambert W Function (hp42s)
(10012020 12:33 PM)Albert Chan Wrote: I don't like testing special cases, if I can avoid it. Nor do I; we all have our preferences. I don't like the fact that this solution depends on RAD mode, for instance ;) I don't immediately see a way out of that. 9 1/X 1/X FP X^2 STO+ ST X SQRT is a bit too long for comfort. Perhaps 9 1/X 1/X FP LN1+X Werner 

10012020, 06:25 PM
Post: #35




RE: Lambert W Function (hp42s)
Using sin(PI) is machine dependent and will take some time.
How about adding +ULP (unit in last place) to the initial value? Code:
Web page updated. http://www.finetune.co.jp/~lyuka/technot...w42s.html 

10012020, 10:25 PM
Post: #36




RE: Lambert W Function (hp42s)
Hi lyuka:
I was wrong to suggest removing x=1/e test. e^W(1/e) take many tries to converge. Because of catastrophically cancelled slope, we got something like this: 0.367879441171442 0.25 0.3051544444752 0.335540374899958 0.351461975177244 0.359608250433447 0.363728172008137 ... I would recommend reverting to previous version (with x=1/e test) Besides, previous version is much easier to understand. 

10022020, 05:44 AM
Post: #37




RE: Lambert W Function (hp42s)
Hi, Albert
Confirmed. And 'eW' on my page is reverted to the previous version 1.1. I was jumping to the conclusion. Thanks. 

10022020, 03:02 PM
Post: #38




RE: Lambert W Function (hp42s)
(09302020 05:28 PM)Albert Chan Wrote: y ← y  (y*ln(y)  x) / (ln(y) + 1) It's rather the y+x that is the culprit  or y*ln(y)x, which is zero as you pointed out. LN(y)+1 for y close to 1/e can be calculated precisely as (LN1P being the LN1+X function) LN(y) + 1 = LN((1/e)*(1+e*(y1/e)) + 1 = 1 + LN1P(e*(y1/e)) + 1 = LN1P(e*(y1/e)) eg. y = 1/e + 1e17 then LN(y) + 1 = 2.71828182845904521e17 LN1P(e*(y1/e)) = 2.718281828459045198415006976699411e17 But the cancellation happens also in y + x, and there's nothing we can do. Or at least, nothing *I* can do ;) Cheers, Werner 

10022020, 05:55 PM
(This post was last modified: 10032020 10:24 AM by Albert Chan.)
Post: #39




RE: Lambert W Function (hp42s)
(10022020 03:02 PM)Werner Wrote: eg. y = 1/e + 1e17 Above amazing accurate results are misleading. You have y1/e = 1e17 (exactly). In other words, y has infinite number of digits ... For rounded 34 digits of y, Free42 will evaluate to the same result. However, both answers only matched half precision. (log1p version does not help here) 1 [EXP] 1e17 [+] // 3.678794411714423315955237701614609e1 = y To get an accurate slope, log(y)+1 = log1p(ε = e*y1) = log1p(ε = e*(y1/e)) But, this shift the problem to get accurate ε, which required more precise 1/e 1e17 // y  1/e 3.255418886896823216549216319830254E35 // (more precise 1/e)  1/e − // 1.000000000000000003255418886896823e17 = y  (more precise 1/e) 1 [EXP] × // 2.718281828459045244209433475626668e17 = ε [LN1+X] // 2.718281828459045207264152980973417e17 = log1p(ε) Mathematica gives 2.718281828459045207264152980973418e17 Quote:But the cancellation happens also in y + x, and there's nothing we can do. Or at least, nothing *I* can do ;) Since y ≈ x, (y + x) is exact, without loss of precision (both x, y are inputs, thus considered exact) 

10022020, 07:29 PM
Post: #40




RE: Lambert W Function (hp42s)
Here is a demo of previous post results, r + err_r = 1/e (108bits accurate)
Code: r, err_r = 0.36787944117144233, 1.2428753672788363E17 >>> from mpmath import * >>> ulp = 2**54 >>> x = r + ulp # 0.36787944117144228 >>> eW(x, g=g0) # simple version, half precision, as expected. 0 0.367879447562281 1 0.367879447022175 2 0.367879445804813 3 0.367879448020601 4 0.367879446543632 mpf('0.36787944701838976') >>> eW(x, g=g1) # log1p version 0 0.367879447562281 1 0.367879447562281 2 0.367879447562281 3 0.367879447562281 4 0.367879447562281 mpf('0.36787944756228136') >>> eW(x, g=g2) # log1p + more precise 1/e 0 0.367879447562281 1 0.367879446846838 2 0.367879446801743 3 0.367879446801563 4 0.367879446801563 mpf('0.36787944680156281') >>> exp(lambertw(x)) # accurate eW mpf('0.36787944680156281') 

« Next Oldest  Next Newest »

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