(28/48/50) Lambert W Function

04012023, 05:36 PM
Post: #21




RE: (28/48/50) Lambert W Function
That sounds very promising. It will take me a while to wrap my head around it though.
As a bonus, we have an accurate Lambert W program for the HP71B. Meanwhile... 

04012023, 05:59 PM
Post: #22




RE: (28/48/50) Lambert W Function
Here is the current version of the Lambert W program for all RPL calculators. This version covers all representable numbers, real or complex, in any branch. In almost all cases the result will be accurate within 2 ULP's, except for a circle of radius .001 around 1/e.
This program replaces the first and last programs in post #1. As before, level 2 should have an integer representing the branch, and level 1 should have the real or complex number z. Code:
Also, a shorter version for the HP 49 and 50 using UNROT, PICK3 and IFTE. The program can be used in approximate or exact mode, but the result will be approximate. Execution time on my 50g ranges from 110 to 490 ms depending on input. Code:


04022023, 11:12 PM
Post: #23




RE: (28/48/50) Lambert W Function
(03262023 06:43 PM)Albert Chan Wrote: Because solving for y = e^W(a) is easier, W code is not recommended. W code is perfect for other branches! (expW code for W branch 0, 1) Numbering the branches of the Lambert W function Except for 2*k*pi*I term, it is basically same formula for W code! lua> f = fn'x: x + I.log(x)  I.log(a)  2*k*pi*I' lua> df = fn'x: 1 + 1/x' Let X = x, A = a log(x)  log(a) = log(X/A) + I*(arg(x)  arg(a)) log(X/A) = 1 + log(X*(r+err)/A) ≈ 1 + log1p((X*err  (A+r) + r*(1+X)) / A) This means we can evaluate f accurately around branch point! f=0 guaranteed we get the correct branch, which is nice. There is no need to use complicated W guesses. lua> a, k = 3+4*I, 5 lua> x = 2*k*pi*I  guess for W_{k}(a) lua> repeat h=f(x)/df(x); x=xh; print(x, I.abs(h)) until x == x+1e6*h (1.815554248884793+30.714634540472343*I) 1.9462907525577031 (1.817005891456528+30.71333413897562*I) 0.0019489253984594729 (1.8170058918466274+30.713334137004896*I) 2.0089622485451688e009 (1.8170058918466274+30.713334137004896*I) 0 >>> from mpmath import * >>> lambertw(3+4j, k=5) mpc(real='1.8170058918466274', imag='30.713334137004896') 

04032023, 07:24 PM
Post: #24




RE: (28/48/50) Lambert W Function
(04022023 11:12 PM)Albert Chan Wrote: f=0 guaranteed we get the correct branch, which is nice. Thanks for the John D. Cook link, I hadn't seen that one before. My program gets exact or nearexact results for all branches, but it is slower for branches other than 0 and 1. I will try your code and see how it compares. 

04032023, 08:47 PM
(This post was last modified: 04032023 09:15 PM by Albert Chan.)
Post: #25




RE: (28/48/50) Lambert W Function
(04032023 07:24 PM)John Keith Wrote: My program gets exact or nearexact results for all branches, but it is slower for branches other than 0 and 1. I was getting stupid. Of course you get nearexact result for other branches. For branches besides (0, 1), a ≈ 1/e is not that important, since W(1/e) ≠ 1 With huge complex parts in the way, there is no catastrophic cancellation issues. (*) The singularity is at only at 0, not 1/e W_{1} is the "hard" branch, with singularities at both places, 0 and 1/e. Plot[{Re[LambertW[1,x]], Im[LambertW[1,x]]}, {x,1,1}] (*) Except for W_{1}: lambertw(z, k) == conj(lambertw(conj(z), k)) >>> lambertw(.36787944+1e99j, k=1) mpc(real='1.0000798057615958', imag='3.4063941215654215e95') >>> conj(lambertw(.367879441e99j, k=1)) mpc(real='1.0000798057615958', imag='3.4063941215654215e95') 

04032023, 10:47 PM
(This post was last modified: 04042023 12:19 AM by Albert Chan.)
Post: #26




RE: (28/48/50) Lambert W Function
(04022023 11:12 PM)Albert Chan Wrote: lua> a, k = 3+4*I, 5 With Newton's rapid convergence, it might be overkill for better W guess ... x = r cis(θ) e^x = e^(r*cos(θ)) cis(r*sin(θ)) x*e^x = r*e^(r*cos(θ)) cis(θ+r*sin(θ)) = a r*e^(r*cos(θ)) = a θ+r*sin(θ) = arg(a) + 2*k*pi = T // for W branch k cos(θ) = ln(a/r)/r sin(θ) = (Tθ)/r We assume k positive > both r and x imag part positive too. lua> a, k = 3+4*I, 5 lua> A, T = I.abs(a), 2*k*pi + I.arg(a) lua> A, T 5 32.34322175389954 lua> c = 0  = cos(θ) lua> r = (Tacos(c)) / sqrt(1c*c)  = T  pi/2 lua> c = log(A/r)/r lua> x = r * I.new(c, sqrt(1c*c)) lua> x  W_{k}(a) guess (1.81718109820366+30.71872424960138*I) Rinse and repeat ... lua> r = (Tacos(c)) / sqrt(1c*c) lua> c = log(A/r)/r lua> x2 = r * I.new(c, sqrt(1c*c)) lua> x2 (1.8170057678766207+30.7133303234809*I) lua> r = (Tacos(c)) / sqrt(1c*c) lua> c = log(A/r)/r lua> x3 = r * I.new(c, sqrt(1c*c)) lua> x3 (1.8170058919343504+30.713334139703402*I) lua> aitken(x, x2, x3) (1.8170058918466334+30.71333413700532*I) We don't even need machine to support complex numbers. We could just iterate for (r, c), improve with Aitken's deltasquared process Assuming K is positve, after (r, c) converged, W_{k}(a) = complex(r*c, r*sqrt(1c*c)) 

04042023, 01:03 AM
(This post was last modified: 04042023 02:26 PM by Albert Chan.)
Post: #27




RE: (28/48/50) Lambert W Function
(04032023 10:47 PM)Albert Chan Wrote: cos(θ) = ln(a/r)/r We can setup solver to get r, then W_{k}(a) Let c = cos(θ), s = sin(θ) s = (Tθ)/r θ = Tr*s = T  r*√(1c²) s ≤ 1 → min(r) = Tθ > Tpi lua> g = fn'r,c: c=log(A/r)/r; acos(c)  (Tr*sqrt(1c*c))' lua> S = require'solver' lua> r = S.secant(g, Tpi, T+pi, nil, true) 29.20162910030975 35.48481440748934 30.766823278237727 30.76703434746101 30.767034374835603 30.767034374835603 lua> c = ln(A/r)/r lua> r * I.new(c, sqrt(1c*c))  = W_{5}(3+4*I) (1.8170058918466274+30.713334137004896*I) I setup function comparing angles. For k ≥ 2, r only have 1 real root. 

04042023, 07:09 PM
Post: #28




RE: (28/48/50) Lambert W Function
Your HP71 program from post #19 is very good. It outperforms my program in both branches for 1/e < x <= .2
I frankly don't feel that FNL from post #20 is worth the extra memory requirement since the post 19 program is always within 3 ULP's. Here is my RPL translation. It is optimized for size rather than readability (lots of "stackrobatics") so I have added comments corresponding to the BASIC statements. As with the other programs, the level 2 argument is k, and the level 1 argument is x. This program is not intended for x > .1, nor for k other than 0 or 1. I'm not sure whether I will integrate it into the main program, which already takes over 600 bytes. For now I'll just leave it asis for anyone who needs higher accuracy near the branch point. Code:


01292024, 11:04 AM
Post: #29




RE: (28/48/50) Lambert W Function
In your post 18,
Hi, Albert. You write 70 Y=(K=1)+K*A/4 80 REPEAT @ H=Y(Y+A)/(LOG(Y)+1) @ Y=YH @ UNTIL Y=Y+H*.0001 What is the meaning of "(K=1)" in line 70 " a) Y=k+k*abs(a)/4 in any case b) or Y=k*abs(a)/4; Y=Y+(1, if k=1; 0, else) c) or Y=? 

01292024, 02:47 PM
(This post was last modified: 01292024 02:58 PM by Gil.)
Post: #30




RE: (28/48/50) Lambert W Function
I tried to translate your code in post 17, Albert.
{ "« x .1 <= IF THEN .367879441171 'R' STO 4.42321595524E13 'R2' STO 'sqrt(2.*R*(x+R+R2))*(k+k+1.)' —>NUM 'Y' STO 'R+Y*sqrt (1.+Y/(3.*R))' —) NUM 'Y' STO DO Y 'Y+x' '(YRR2)/R' —>NUM DUP TYPE 1. == IF THEN 1. + LN ELSE LNP1 END /  —>NUM 'H' STO 'YH' —>NUM 'Y' STO Y 'Y+H*.0001' —>NUM == UNTIL END ELSE 'k+k+1.' '(k+k+1.)*x/4.' + —>NUM 'Y' STO DO 'Y(Y+x)/(LN(Y)+1.)' —>NUM 'H' STO 'YH' —>NUM 'Y' STO UNTIL Y 'Y+H*.0001' —>NUM == END END x Y / »x =.36787944118 W(k=>; (.999999999985,6.87304618656E6) Result ± ok for imaginary part. Your translated code for more accuracy, post 19, same input values « x .1 IF THEN .367879441171 'R' STO 4.42321595524E13 'R2' STO '(x+R+R2)/R' —>NUM 'H' STO 'sqrt(2.*H)*(2.*k+1.)' —>NUM 'X' STO 'X*sqrt(1.+X/3.)' —>NUM 'X' STO 'R+R*X' —>NUM 'Y' STO DO X DUP TYPE 1. == IF THEN 1. + LN ELSE LNP1 END 'Z' STO 'X(XZ+H)/Z', —>NUM 'Z' STO X Z  'X' STO X 'X+Z*.000001' —>NUM == UNTIL END X DUP TYPE 1. == IF THEN 1. + LN ELSE LNP1 END 1.  ELSE 'k+k+1.' '(k+k+1.)*x/4.' + —>NUM 'Y' STO DO 'Y(Y+x)/(LN(Y)+1.)' —>NUM 'H' STO 'YH' —>NUM 'Y' STO UNTIL Y 'Y+H*.0001' —>NUM == END x Y / END » No solution is found: X is never equal to 'X+Z*.000001'. Why this difference of "behaviour"? And why, in code post 29, did you define Y=R+R*X, that seems never to be used again? 

01292024, 03:51 PM
(This post was last modified: 01292024 05:20 PM by Albert Chan.)
Post: #31




RE: (28/48/50) Lambert W Function
Hi, Gil
Post 18 HP71B code is short, I'll just copy/pasted in full Code: 10 INPUT "a, k? ";A,K @ K=K+K+1 Line 10: K=K+K+1 map W branch (0,1) into sign ±1 Line 40: this simplified lyuka's formula 2 branches ±√ term But lyuka's formula is not suitable if A is too far from 1/e Example, if √ term keep growing, eventually we would get e^W guess = 0. But, this does not match exactly at A = 0, which make W_{1} guess very bad. For that, I just use simple guess (note: K = ±1) Line 70: Y=(K=1)+K*A/4 K=1 > Y = A/4 // > X = ln(A/4), roughly ln(A) if A → 0 K=+1 > Y = 1 + A/4 // > X = log1p(A/4), roughly ln(A) if A → INF There was another way to get W, based from how lyuka's formula work? Then, solve the same thing using accurate log(1+x)x Turns out, this super accuracy is not needed. I made a version with FNL(x) = accurate log(1+x)x, and it make little difference. (post #20) (03312023 10:07 PM)Albert Chan Wrote: When we are close to branch point, x is tiny. Even rough x estimate suffice. Gil Wrote:No solution is found: X is never equal to 'X+Z*.000001'. Math assumed real numbers to begin with. It may work with complex numbers, but you still need adjustments. 

01292024, 06:46 PM
Post: #32




RE: (28/48/50) Lambert W Function
Thanks.
And why, in code post 29, did you define Y=R+R*X, that seems never to be used again? 

01292024, 07:30 PM
Post: #33




RE: (28/48/50) Lambert W Function
(01292024 06:46 PM)Gil Wrote: And why, in code post 19, did you define Y=R+R*X, that seems never to be used again? (R+R*X) is just another method to calculate Y = e^W(A) It is only for compare against Newton method Y = (Y+A)/LOGP1((YRR2)/R) I was trying to decide which way does the job better. (12302022 11:27 PM)Albert Chan Wrote: When I try to understand lyuka's eW formula (previous post), I found a better one. I am skipping details of how to estimate Z from H. Once we have Z, e^W(a ≈ 1/e) problem is solved. e^W(a) = y = r + z = r + Z/e = r + r*Z // = R+R*X of HP71B code 

01292024, 09:50 PM
Post: #34




RE: (28/48/50) Lambert W Function
What hawk eyes, Albert! Not only your brain works marvels, Albert, but your eyes too.
Sorry for having mixed up the post number! More seriously: for k=1, and a=.36787944118 (>1/e, about .367879441171), how do I know that formulae post 19, more accurate, will fail (at least with the code I put in my calculator) and that formulae post 18, a wee bit less accurate, will work? In other words: should I say, for "a" near 1/e, but just smaller — as in the example above —, I should always play safety and only use formulae in post 18? Second question k mappin — new k = ± k. Does it mean that the derived formulae in posts 18 and 19 only work for k, initially = 0 or 1? If yes, are there other pitfalls to look at for initial k≠{0, +1} when "a" is near 0 or 1/e? Regards 

01302024, 12:18 AM
Post: #35




RE: (28/48/50) Lambert W Function
(01292024 09:50 PM)Gil Wrote: for k=1, and a=.36787944118 (>1/e, about .367879441171), how do I know that formulae post 19, more accurate, will fail (at least with the code I put in my calculator) and that formulae post 18, a wee bit less accurate, will work? They are equally accurate, if coded correctly. They are really the same e^W formula, except we let Y = R+R*X It is just scale and offset, convergence rate are exactly the same. There is no wrong answer. Pick your preferred choice. My preference is Y = (Y+A)/LOGP1((YRR2)/R), because convergence is easily tested. It is also easy to explain. Denominator is simply more accurate (LN(Y)+1) The other version, X = (XL+H)/L, where L = LOGP1(X), is harder to know when to stop. Because of inaccurate XL, X itself never really converged, only (1+X) = Y/R does. OTTH, X formula is simpler than Y's (also, (R, R2) not used at all!) With limited domain (A ≈ 1/e), fixed loops is simple to implement. Quote:If yes, are there other pitfalls to look at for initial k≠{0, +1} when "a" is near 0 or 1/e? Yes! Only used this for k=0, or ±1 with small_imag part (k and im(a) of opposite sign) This is because y = e^x lost branch information. Good guess is crucial to get right branch! e^(2*k*pi*I) = cis(2*k*pi) = 1 // k is gone What happen if k = ±1 and, k and im(a) have same sign? x + ln(x) = ln_{k}(a) im(x) + arg(x) = (arg(a) + 2*k*pi) = [±2*pi .. ±3*pi] im(x) = [±2*pi .. ±3*pi]  [0 .. ± pi] = [± pi .. ±3*pi] Away from real line, no singularity around 1/e, easy to solve. 

01302024, 12:33 AM
Post: #36




RE: (28/48/50) Lambert W Function
Thanks.
Strange anyway, a = .36787944118, that in my calculator the until equality is never reached with method exposed in post 11.—> endless loops. 

01302024, 01:09 AM
Post: #37




RE: (28/48/50) Lambert W Function
(01302024 12:18 AM)Albert Chan Wrote: With limited domain (A ≈ 1/e), fixed loops is simple to implement. For my Python code (~ 15 digits precision), it never get over 5 loops. Perhaps you can use fixed Newton loops instead of testing convergence. Also, you may use this to implement complex log1p(x) y=1+x > log1p(x) ≈ ln(y)  (y1x)/y https://www.hpmuseum.org/forum/thread55...#pid111682 

01302024, 12:04 PM
(This post was last modified: 01302024 12:57 PM by Gil.)
Post: #38




RE: (28/48/50) Lambert W Function
Hi, Albert.
With your function definition of post 20, 1) could I modify as follows? 200 DEF FNL(X) ! « —> X «210 IF ABS(X)>=.4 THEN X=X/(SQRT(1+X)+1) @ L=FNL(X)*2X*X ELSE 220 X2=X/(X+2) @ X4=X2*X2 230 X4=X4*(5005X4*(5082X4*969))/(15015X4*(24255X4*(11025X4*1225))) 240 L=X2*(X4+X4X) 250 END » » for my HP50G calculator program. 2) Could I, consequently, simplify "your" line 42 as follows? 42 DO FNL(X) @Z=X(HL)/(X+L) @ X=XZ @ UNTIL X=X+Z*.0001 END 44 PRINT "e^W, W =";R+R*X;LOGP1(X)1 3) And I should start the program with line 40 defined as before >40 H=(A+R+R2)/R @ X=SQRT(2*H)*K @ X=X*SQRT(1+X/3) @ Y=R+R*X Correct? In fact, it seems that L cannot be equal to @ L=FNL(X)*2X*X Or the defined function should be rather 200 DEF FNL(X) ! « —> X «210 IF ABS(X)>=.4 THEN X=X/(SQRT(1+X)+1) @ L=X*2X*X ELSE 220 X2=X/(X+2) @ X4=X2*X2 230 X4=X4*(5005X4*(5082X4*969))/(15015X4*(24255X4*(11025X4*1225))) 240 L=X2*(X4+X4X) 250 END Thanks again, Albert, for your possible and always precious help. Gil 

01302024, 01:47 PM
Post: #39




RE: (28/48/50) Lambert W Function
Hi, Gil
FNL(x) = accurate ln(1+x)x, was designed for real number. It may need work to extend for complex numbers. Also, it is defined recursively, until x is small, not a IF THEN ELSE structure. see thread: Accurate x  log(1+x) (for consisteny, I switched order, and do log1p_sub(x)) Anyway, accurate FNL(x) is not needed here. Here is a direct replacement snippet for my Python W code I am using formula to get log1p, instead of calling builtin log1p Quote:y=1+x > log1p(x) ≈ ln(y)  (y1x)/y Because of inaccurate xlog1p(x), x does not converge, but (1+x) does. (again, you may forget about checking convergence, and just loop a few times) Code: if abs(h)<.25 and (k==0 or small_imag): With identical termination condition, results compared well to W original version. Code: if abs(h)<.25 and (k==0 or small_imag): 

01302024, 02:52 PM
Post: #40




RE: (28/48/50) Lambert W Function
I tried and changed the code fo my function FNL relative to ln(1+x)x
FNL= « WHILE x .4 >= REPEAT 'x/(sqrt(1+x)+1)' >NUM 'x' STO x FNL 2 * x SQ  'fnl' STO END 'x/(x+2)' >NUM 'X2' STO X2 DUP * 'X4' STO 'X4*(5005X4*(5082X4*969))/(15015X4*(24255X4*(11025X4*1225)))' >NUM 'X4' STO 'X2*(X4+X4x)' >NUM 'fnl' STO » But, for initial x=4, I get a final fnl = 4.77474540488E2 that is incorrect, of course. 

« Next Oldest  Next Newest »

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