(28/48/50) Lambert W Function

03202023, 08:43 PM
(This post was last modified: 06282023 06:23 PM by John Keith.)
Post: #1




(28/48/50) Lambert W Function
NOTE: The first and last programs have been obsoleted by the program in post #22, which should be used in most cases. That program gives accurate results over the entire complex plane for any branch, positive or negative, but accuracy is limited for values very close to the branch point 1/e.
For very accurate results near the branch point in branches 0 and 1, use the program in post #28, which is a translation of Albert Chan's HP71 program in post # 19. The following text refers to the (mostly obsolete) programs in this post only.  This program computes the Lambert W function W(z) for branches 0 and 1. For branch 1, the program covers the realvalued range 1/e < z < 0. For branch 0, the program covers the entire complex plane. For most values, the results are accurate to within 2 ULPs. For values very close to 1/e, the results are approximate. The program expects the branch (0 or 1) on level 2 and z on level 1. Some notes and background information: The gold standard for Lambert W programs on the HP 50g is Egan Ford's SPECIAL.LIB. It is very accurate for almost any input value, and can return results for any branch. It does have some problems, however. It is an HPGCC program and only works on the 49g+ and 50g. It requires two large programs on the calculators SD card. Also, it is rather slow. I believe that the slowness is caused by the need to fetch the executable from the SD card, as one would expect the HPGCC code to be quite fast. The program described here is comparatively limited since no branches other than 0 and 1 are covered, and then only the realvalued range of branch 1. However, it is about 5 times as fast as SPECIAL.LIB and about 9 times as fast as programs using the calculator's builtin solvers. It also runs on all RPL calculators. If anyone here has an idea of how to get an accurate estimate for other branches, or to improve accuracy around the branch point at 1/e, please chime in! The code that provides the initial estimate for branch 0 comes from page 8 of Iacono and Boyd. The estimation code for branch 1 comes from formula 27 on page 12 of https://www.sciencedirect.com/science/ar...via%3Dihub. The recurrence used in the main loop occurs in both papers and is mentioned in Wikipedia. Many additional references here. Code:
Also, a streamlined version for the 49 and 50, 25 bytes smaller: Code:
Now, here are a couple of experimental programs for those interested in further exploration. The first one is a variation of the recurrence from the main program using the LongFloat Library. The program expects z on level 2 and an estimate (such as the result of the main program) on level 1. It is rather slow, about 22 seconds on my 50g but reasonable on emulators. The value of c assumes the default value of DIGITS which is 50. Removing the R\<\>F at the end will leave a LongFloat number (type 27) on the stack. Code:
The last program implements the Halley recurrence shown at the Wikipedia link and in the Lóczi paper. The reason for including it is that the recurrence used in the main program works poorly for branches other than branch 0 and the realvalued part of branch 1. It tends to "drag" the values back to branch 0 even if given an accurate estimate for the value in another branch. The Halley recurrence converges rapidly even if the estimate is not very close, although it can be fooled if the estimate is vague. As above, the program expects z on level 2 and an estimate (usersupplied) on level 1. For HP28 and 48, change PICK3 to 3 PICK. Code:


03202023, 10:51 PM
Post: #2




RE: (28/48/50) Lambert W Function
The issue is how to accurately calculate Newton's iteration, around branch point, 1/e
I prefer solving for y = e^W(a). Its Newton's method is very simple. x*e^x = ln(y)*y = a f = y*ln(y)  a y = y  f/f' = y  (y*ln(y)  a) / (ln(y)+1) y = (y+a) / (ln(y)+1) Around a ≈ 1/e, (y, a) are considered input, thus "exact". They are about the same size, with opposite sign. By Sterbenz lemma, (y+a) is exact. e^W(1/e) = 1/e // a ≈ 1/e > y ≈ +1/e (12302022 11:48 PM)Albert Chan Wrote: ln(y)+1 = ln(e*y) = log1p(e*(y1/e)) For accurate denominator, we need "doubled" precision for 1/e For 12 digits, 1/e ≈ 0.367879441171 + 0.442321595524E12 To reduce iterations, we also need good starting estimate. Lua e^W code (0, 1 branch): https://www.hpmuseum.org/forum/thread19...#pid167919 

03212023, 06:15 AM
Post: #3




RE: (28/48/50) Lambert W Function
This programme works for positive reals  I'm interested in the comparison of accuracy:
Code: :: 

03212023, 01:53 PM
Post: #4




RE: (28/48/50) Lambert W Function
(03202023 10:51 PM)Albert Chan Wrote: For accurate denominator, we need "doubled" precision for 1/e That's pretty much what I feared. I imagine that the best compromise would be to use extended precision reals, but SysRPL is beyond my skill set. The estimation code I used for branch 0 is amazingly accurate over the entire numerical range representable by Saturnbased HP's. Also, the recurrence that I used converges faster than Newton's method while not being much larger in code size. Newton's method may be somewhat more accurate near the branch point but it wouldn't be much of an improvement without doubleprecision numbers. Of more interest to me is computing W(z) for branches other than 0 and 1. I have not been able to find any methods for doing so. Any insight in this area would be greatly appreciated. 

03212023, 05:15 PM
Post: #5




RE: (28/48/50) Lambert W Function
(03212023 01:53 PM)John Keith Wrote: I imagine that the best compromise would be to use extended precision reals, but SysRPL is beyond my skill set. Just to be clear, we don't need extended precison numbers. My Lua code only use the same 64bits float throughout. All we need is more precise constant 1/e, so that (y1/e) does not suffer cancellation errors. Quote:The estimation code I used for branch 0 is amazingly accurate over the entire numerical range representable by Saturnbased HP's. Also, the recurrence that I used converges faster than Newton's method while not being much larger in code size. I am unfamilar with RPL, can you post actual formula for W(a) guess, and the iteration formula? Convergence can be affected greatly by what the iteration formula look like. Newton's method is not necessarily much slower than higher order method. To give an example, lets do W(a = 1e10), with guess a0 = log1p(a) ≈ 23.0259 Note: nestlist is my defined function, not part of XCas XCas> a := 1e10; a0 := log1p(a) XCas> f := x*e^x  a XCas> nestlist(unapply(xf/f',x), a0, 3) → [23.0259, 22.1091, 21.2606, 20.568] XCas> nestlist(unapply(xf/(f'f''/2*f/f'),x), a0, 3) → [23.0259, 21.2714, 20.1788, 20.029] This is a bad setup, with huge derivatives, f' = (x+1)*e^x, f'' = (x+2)*e^x If a is bigger , bad setup is going to take longer to converge. Let's try a better setup. XCas> f := x  ln(a/x) XCas> nestlist(unapply(xf/f',x), a0, 3) → [23.0259, 20.0198, 20.0287, 20.0287] XCas> nestlist(unapply(xf/(f'f''/2*f/f'),x), a0, 3) → [23.0259, 20.0279, 20.0287, 20.0287] Here, there is not much need for complicated halley's method. Newton is almost as good. If a is bigger , better setup is going to converge faster! Perhaps in 1 iteration! Quote:Of more interest to me is computing W(z) for branches other than 0 and 1. I have not been able to find any methods for doing so. Any insight in this area would be greatly appreciated. Python mpmath lambertw(z, k=0) support all branches. https://mpmath.org/doc/0.19/functions/po...wfunction https://github.com/mpmath/mpmath/blob/ma...ns.py#L464 

03222023, 07:30 PM
Post: #6




RE: (28/48/50) Lambert W Function
(03212023 05:15 PM)Albert Chan Wrote: I am unfamilar with RPL, can you post actual formula for W(a) guess, and the iteration formula? Many thanks for the mpmath github reference, that solved my branch problem! Updated program will follow soon. Here are the formulas that you requested: Formula for global approximation of W(z) for branch 0. From "New approximations to the principal realvalued branch of the Lambert Wfunction", by Iacono and Boyd. 1 + a*ln((1 + b*y)/(1 + c*ln(1+y)) where y = sqrt(1+e*x) c = (e^(1/a)  1  sqrt(2)/a)/(1  e^(1/a)*ln(2)) b = sqrt(2)/a + c With 'a' determined empirically to be 2.036, final formula is: 1 + 2.036*ln((1 + 1.14956131*y)/(1 + 0.4549574*ln(1 + y)) Recurrence from same paper, also from "Guaranteed and highprecision evaluation of the Lambert W function" by Lajos Lóczi. B'(x) = (B(x)/(1 + B(x))*(1 + ln(x/B(x)) where B(x) is current estimate B'(x) is new estimate 

03222023, 08:55 PM
(This post was last modified: 03222023 08:56 PM by John Keith.)
Post: #7




RE: (28/48/50) Lambert W Function
(03212023 06:15 AM)Gerald H Wrote: This programme works for positive reals  I'm interested in the comparison of accuracy: They are roughly comparable, unless x is very large. A few examples: Code:
Exact values were computed with Mathematica version 12.2. As you can see, both programs lose accuracy for values very close to 1/e. Your program is very accurate over most of the real range but less accurate at MAXR. This is a known problem with Newton's method. 

03232023, 12:16 AM
(This post was last modified: 03232023 12:17 AM by Albert Chan.)
Post: #8




RE: (28/48/50) Lambert W Function
(03222023 07:30 PM)John Keith Wrote: y = sqrt(1+e*x) There is a flaw with Iacono and Boyd W guess formula. If x is tiny enough, below 6e10, it could return a negative guess. This cause Newton formula to take a log of negative number, which crash it. 

03232023, 02:56 AM
(This post was last modified: 04042023 01:41 PM by Albert Chan.)
Post: #9




RE: (28/48/50) Lambert W Function
Here is my implementation of accurate x = W(a), both branches.
x * e^x = a ln(x) + x = ln(a) → f = x + ln(x/a) = 0 x  f/f' = x  (x + ln(x/a)) * x/(1+x) = (1  ln(x/a)) * x/(1+x) W(1/e) = 1 > a ≈ 1/e, x ≈ 1 > (1+x) is tiny Let 1/e = r + err 1  ln(x/a) =  ln(x*(r+err)/a) =  log1p((x*err(a+r)+r*(1+x))/a) x + ln(x/a) = (x+1)  (1  ln(x/a)) Terms inside log1p sorted in increasing size, r*(1+x) the biggest. With this, and some minor changes to my expW function, I get this: Code: function W(a, x, verbal) lua> W(0.367, 0, true) 0.9323991847479226 0.9323991847479282 lua> W(0.367, 1, true) 1.0707918867680617 1.0707918867680521 lua> W(1e10, 0, true) 18.111645253927524 20.02372402639066 20.028685384073526 20.02868541330495 20.028685413304952 Update 1: added special case W_{0}(0) = 0, W_{1}(0) = ∞ (should it be NaN?) Update 2: removed guard "if h == err then return 1 end", see next post. Update 3: improve W guess with a Newton's step, based from y = e^W y = (y+a) / (log(y)+1) → x = a/y = a/(y+a) * (log(y)+1) Update 4: change behavior from "W branch 1 if x" to "x default to branch 0" Assumed input, x = nil/false/0/1 only: x = nil/false/0 get branch 0; x = 1 get branch 1 

03232023, 05:53 PM
Post: #10




RE: (28/48/50) Lambert W Function
(03232023 12:16 AM)Albert Chan Wrote: There is a flaw with Iacono and Boyd W guess formula. That explains the problem that I saw when z was < 10^8. Lines 3 and 7 of my first program are a patch around the flaw. The program doesn't crash without the patch but it returns a complex result in a different branch. There is one line in your program that concerns me: if h == err then return 1 end It seems that the program would return 1 if a is the closest machine representation of 1/e, which would be .367879441171 on the HP calculators. Instead, we would want the program to return .999998449288 in branch 0, or 1.00000155071 in branch 1, which are the correct 12digit representations of W(.367879441171). My apologies if I am reading your program incorrectly. Also, there is another line whose meaning is unclear to me: x = sqrt(2*h/r) * (x and 1 or 1) In RPL this statement always returns 1. 

03232023, 07:30 PM
(This post was last modified: 03232023 07:42 PM by Albert Chan.)
Post: #11




RE: (28/48/50) Lambert W Function
(03232023 05:53 PM)John Keith Wrote: There is one line in your program that concerns me: You are right. Guard should be if (h == 0) ... but that would never happen. Guard is not needed. I'll remove it. Quote:x = sqrt(2*h/r) * (x and 1 or 1) Last term is Lua idiom for ternary operators It says that pick 1 if x evaluated true (neither nil nor false), 1 otherwise. Signs are the result of solving W guess quadratic equation, 2 roots for W 2 branches 

03262023, 06:43 PM
Post: #12




RE: (28/48/50) Lambert W Function
(03232023 02:56 AM)Albert Chan Wrote: Here is my implementation of accurate x = W(a), both branches. I just realized this is equivalent to solving for y = e^W(a) With "same" guess, x*y=a, Newton convergence rate are identical. x = (1  ln(x/a)) * x/(1+x) a/y = (1 + ln(y)) * 1/(1+y/a) y = (y+a) / (1 + ln(y)) lua> expW(1e99, nil, true) 2.5e+098 4 5.492824331831047e+096 182.05570387623288 4.493790313227532e+096 222.52929716290643 4.493356750520384e+096 222.55076895111614 4.493356750426822e+096 222.55076895575016 4.493356750426821e+096 222.55076895575021 lua> W(1e99, nil, true) 182.05570387623249 222.52929716290646 222.55076895111617 222.55076895575016 222.5507689557502 Note: W guess used 1 Newton step. Both code converged the same way. Because solving for y = e^W(a) is easier, W code is not recommended. Bonus: with y, there are 2 ways to recover W(a) = a/y or log(y) Depends on size of y (use log if y is big), we can make W(a) estimate better. Example, above last 2 y's differ by 1 ULP, but the log's are the same. (Not exactly. Here, y 1 ULP error translated to 0.007 ULP error in x) lua> log(4.493356750426822e+096) 222.5507689557502 lua> log(4.493356750426821e+096) 222.5507689557502 

03272023, 04:45 PM
Post: #13




RE: (28/48/50) Lambert W Function
Hi Albert,
I am still having problems with your Lua program. This line, h = x/(1+x) * s(x) returns zero or a very small number (usually 10^12 or less) and thus no correction occurs in the line x = x  h. This may be caused by the difference between the 64bit binary floats used by Lua and the 12digit BCD floats used by HP calculators. It is also possible that I made an error in translating the program although I have checked carefully. I was careful to sum the values in this line s = function(x) return (x+1) + log1p((x*err(a+r)+r*(1+x))/a) in the proper order as suggested in your post. I have also seen several mentions of using Taylor series rather than iteration for values close to 1/e. This may also be worth exploring. 

03272023, 08:57 PM
(This post was last modified: 03272023 11:21 PM by Albert Chan.)
Post: #14




RE: (28/48/50) Lambert W Function
(03272023 04:45 PM)John Keith Wrote: This line, h = x/(1+x) * s(x) returns zero or a very small number (usually 10^12 or less) Hi, John Keith I assumed you meant correction around branchpoint is small. That's normal. My W guess formula is pretty accurate for a ≈ 1/e, both branches. Solve s(x) = x + ln(x/a) = 0 (or its messy version) is probably a bad idea. My previous post suggested it is better to solve for y = e^W(a), shown below. lua> a = 0.3678  calculate W(a), 0 branch, by "hand" lua> r, err = 0.36787944117144233, 1.2428753672788363E17 lua> h = a + r + err lua> x = sqrt(2*h/r)  rough guess for e*h = (1+x)*log1p(x)  x lua> x = x * sqrt(1+x/3)  better guess lua> y = r + r*x  = e^W(a), if x is correct lua> y, a/y, log(y) 0.3755511062373703 0.9793607152032429 0.9793607152084147 W(a) 2 ways matched does not imply they are correct, only that we are close. Interestingly, true W(a) may not be a "mean" of the two ways. lua> y = (y+a) / log1p((yrerr)/r)  Newton's step lua> y, a/y, log(y) 0.37555110633147754 0.9793607149578304 0.9793607149578304 y's had doubled in accuracy, to full precision = e^W(a) Quote:I have also seen several mentions of using Taylor series rather than iteration for values close to 1/e. XCas> H := (1+x)*ln(1+x)  x XCas> x0 := sqrt(2*H) /* rough guess for above formula */ XCas> r0 := revert(series(x0,0,8)) \(\displaystyle x +\frac{x^2}{6} \frac{x^3}{72} +\frac{x^4}{270} \frac{23 x^5}{17280} +\frac{19 x^6}{34020} \frac{11237 x^7}{43545600} +\frac{13 x^8}{102060}\) XCas> rs := series((r0/x)^2) /* correction, inside square root */ \(\displaystyle 1+\frac{x}{3} +\frac{x^3}{360} \frac{x^4}{810} +\frac{23 x^5}{40320} + O\!\left(x^6\right)\) This explained why my W guess formula is good (no x^2 terms, others tiny, with opposite sign) Let's numerically check if series is good, with pade approximation matching to x^4 (replace 57 by 4032/71 ≈ 56.8, pade matched to x^5, but make no difference here) lua> x = sqrt(2*h/r)  rough guess for e*h = (1+x)*log1p(x)  x lua> x = x * sqrt(1+x/3 + x^3/(360+160*x*(1x/57))) lua> y = r + r*x  = e^W(a), if x is correct lua> y, a/y, log(y) 0.3755511063314775 0.9793607149578305 0.9793607149578305 For comparison, my W code use simple guess (no bolded mess) + 1 Newton step. Newton's step is very cheap; we had to convert x to W(a) guess anyway. W(a) = a/y = a * (ln(y)+1)/(y+a), with y = r+r*x, we have: Code: x = a * log1p(x)/(r*x+h)  Newton step lua> W(0.3678, nil, true)  h=0 on the first try. 0.9793607149578305 0.9793607149578305 

03272023, 11:30 PM
(This post was last modified: 03302023 03:58 PM by Albert Chan.)
Post: #15




RE: (28/48/50) Lambert W Function
(03272023 08:57 PM)Albert Chan Wrote: lua> h = a + r + err Technically, it should be y = (1+x)/e = (1+x)*(r+err) Numerically, err < 1/2 ULP of r, thus y = r + r*x e*h = (e*y) * log(e*y)  (e*y1) h = y * (log(y)+1)  y + 1/e h  1/e = y*log(y) = a A novel way is to directly solve for x, using log1p_sub(x), then convert to W(a) or e^W(a) Except for getting accurate h, err = 1/e  float(1/e) is not used anymore. Let H=e*h, L = log(1+x)  x, accurately computed f = (1+x)*log(1+x)  x  H f' = log(1+x) + (1+x)/(1+x)  1 = log(1+x) x  f/f' = (HL) / (x+L) Note that numerator is really an addition, with both H, L positive. lua> x  better guess, from above quote 0.020853747742736108 lua> H = h/r lua> L = log1p_sub(x) lua> x = (HL) / (x+L)  newton step lua> x  converged to full precision 0.020853747998545925 lua> r+r*x, log1p(x)1  = e^W(a), W(a) 0.3755511063314775 0.9793607149578305 

03312023, 04:06 PM
Post: #16




RE: (28/48/50) Lambert W Function
(03272023 08:57 PM)Albert Chan Wrote: lua> x = sqrt(2*h/r)  rough guess for e*h = (1+x)*log1p(x)  x I was curious why the 2 W's do not bracket true W(a) Doing a bit of error analysis, everything becomes clear. Ignoring error of division, we have: relerr(x = a/y) ≈  relerr(y) ln(y*(1+ε)) ≈ ln(y) + ε ln(y*(1+ε)) / ln(y) ≈ 1 + ε/ln(y) relerr(ln(y)) ≈ relerr(y) / ln(y) 2 W's don't bracket true W(a) is because ln(y) ≈ x is negative. 2 W's have errors of the same sign. relerr(x) + x*relerr(ln(y)) ≈ relerr(y) + relerr(y) ≈ 0 lua> x, lny = 0.9793607152032429, 0.9793607152084147 lua> (x + x*lny) / (1+x) 0.9793607149578298 Simply extrapolate for 0 error, we get good W(0.3678) estimate. Actually, extrapolate for 0 error is equivalent to a Newton's step! 

03312023, 06:15 PM
Post: #17




RE: (28/48/50) Lambert W Function
(03312023 04:06 PM)Albert Chan Wrote: lua> x, lny = 0.9793607152032429, 0.9793607152084147 Thanks for your help and insights, this is all quite interesting if a bit beyond my level of knowledge. Unfortunately, I am still not able to achieve that level of accuracy with the 12digit precision of HP calculators. For an input of .3678, your program for e^W(x) returns .375551106194 and .979360715069 Applying your correction above returns .979360714941 which is a noticeable improvement but still off by 16 ULP's. My program from post #1 returns .979360714903 which is a bit worse (off by 54 ULP's) but with no correction applied. What I am hoping for is a way to achieve 11digit accuracy for values very close to 1/e, such as .36787944117. I still don't think this will be possible without extended precision but I am hoping to be pleasantly surprised. I am close to completing a new version of my program which covers the entire complex plane in all branches. It seems to be accurate to within 2 ULP's except for a circle of radius ~.001 around 1/e. It is only that area that i am concerned about. 

03312023, 07:10 PM
Post: #18




RE: (28/48/50) Lambert W Function
(03312023 06:15 PM)John Keith Wrote: Unfortunately, I am still not able to achieve that level of accuracy with the 12digit precision of HP calculators. My Lua expW code translated to HP71B. It seems to work fine. Code: 10 INPUT "a, k? ";A,K @ K=K+K+1 >run a, k? .3678, 0 e^W, W = .37555110633 .979360714962 >run a, k? .3678, 1 e^W, W = .360260737301 1.02092723941 >run a, k? 1e6, 0 e^W, W = .999998999999 .000001000001 >run a, k? 1e6, 1 e^W, W = 6.01449171278E8 16.6265089014 

03312023, 10:07 PM
Post: #19




RE: (28/48/50) Lambert W Function
(03272023 11:30 PM)Albert Chan Wrote: Let H=e*h, L = log(1+x)  x, accurately computed When we are close to branch point, x is tiny. Even rough x estimate suffice. When we are far away, log1p(x)x does not suffer catastrophic cancellation. We may not need accurate log1p(x)x after all. Here, I added code to roughly solve above f=0 for x, then convert to W's Old code with iterations of y remained (2nd e^W, W), for comparison. >40 H=(A+R+R2)/R @ X=SQRT(2*H)*K @ X=X*SQRT(1+X/3) @ Y=R+R*X >42 REPEAT @ Z=LOGP1(X) @ Z=X(XZ+H)/Z @ X=XZ @ UNTIL X=X+Z*.000001 >44 PRINT "e^W, W =";R+R*X;LOGP1(X)1 >run a, k? .367879441171,0 e^W, W = .367880011645 .999998449291 ! true W = .999998449288 e^W, W = .367880011646 .999998449291 >run a, k? .36787944117,0 e^W, W = .367880471317 .999997199776 ! true W = .999997199775 e^W, W = .367880471317 .999997199778 >run a, k? .3678794411,0 e^W, W = .367886691321 .999980292244 ! true W = .999980292245 e^W, W = .36788669132 .999980292247 >run a, k? .367879441,0 e^W, W = .367890672444 .999969470701 ! true W = .999969470701 e^W, W = .367890672444 .999969470702 ... >run a, k? .3678,0 e^W, W = .375551106332 .979360714956 ! true W = .979360714958 e^W, W = .37555110633 .979360714962 

04012023, 12:44 PM
(This post was last modified: 04012023 01:16 PM by Albert Chan.)
Post: #20




RE: (28/48/50) Lambert W Function
(03312023 10:07 PM)Albert Chan Wrote: We may not need accurate log1p(x)x after all. With accurate log1p(x)x, we can get W around branch point to within 1 ULP or less. With below FNL(X), previous post examples all match W exactly, no ULP errors. 200 DEF FNL(X) ! = ln(1+X)  X, but more accurate 210 IF ABS(X)>=.4 THEN X=X/(SQRT(1+X)+1) @ FNL=FNL(X)*2X*X @ GOTO 250 220 X2=X/(X+2) @ X4=X2*X2 230 X4=X4*(5005X4*(5082X4*969))/(15015X4*(24255X4*(11025X4*1225))) 240 FNL=X2*(X4+X4X) 250 END DEF line 230: X4 = pade approx. of atanh(X2)/X2  1 ≈ ln(1+X)/(2*X2)  1 line 240: FNL = X2*(X4+X4X) ≈ ln(1+X)  2*X2  X*X2 = ln(1+X)  X 42 REPEAT @ L=FNL(X) @ Z=X(HL)/(X+L) @ X=XZ @ UNTIL X=X+Z*.0001 44 PRINT "e^W, W =";R+R*X;LOGP1(X)1 Away from branch point is a different story. With X getting big, FNL(X) won't help much. 

« Next Oldest  Next Newest »

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