Lambert W Function (hp-42s)
09-28-2020, 04:06 PM
Post: #21
 lyuka Junior Member Posts: 48 Joined: May 2014
RE: Lambert W Function (hp-42s)
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,
So I wrote an e^W calculation program for 42S, and additionally in C.

http://www.finetune.co.jp/~lyuka/technot...w-42s.html

It's nice to be able to handle complex numbers easily with 42S,
but if you try to find e^W with Newton-Raphson 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.
09-28-2020, 10:01 PM
Post: #22
 Albert Chan Senior Member Posts: 2,103 Joined: Jul 2018
RE: Lambert W Function (hp-42s)
(09-28-2020 04:06 PM)lyuka Wrote:  It's nice to be able to handle complex numbers easily with 42S,
but if you try to find e^W with Newton-Raphson method, it will fail very close to -1/e.

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.54426772287e-11j)
09-29-2020, 09:21 AM (This post was last modified: 09-29-2020 09:30 AM by lyuka.)
Post: #23
 lyuka Junior Member Posts: 48 Joined: May 2014
RE: Lambert W Function (hp-42s)
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 + 1E-11).

Using Free42 with the example value of x = -1 / e + 1E-11 ~= -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.
09-29-2020, 07:07 PM (This post was last modified: 09-29-2020 08:17 PM by Albert Chan.)
Post: #24
 Albert Chan Senior Member Posts: 2,103 Joined: Jul 2018
RE: Lambert W Function (hp-42s)
(09-28-2020 04:06 PM)lyuka Wrote:  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.

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 divide-by-zero.
You might want to adjust the formula a tiny bit ...
09-29-2020, 11:17 PM
Post: #25
 lyuka Junior Member Posts: 48 Joined: May 2014
RE: Lambert W Function (hp-42s)
That case of x = -1 / e has already been handled by my code shown in the link in the previous post.
09-30-2020, 07:08 AM (This post was last modified: 09-30-2020 07:09 AM by Werner.)
Post: #26 Werner Senior Member Posts: 748 Joined: Dec 2013
RE: Lambert W Function (hp-42s)
(09-28-2020 04:06 PM)lyuka Wrote:  http://www.finetune.co.jp/~lyuka/technot...w-42s.html
In that program:

Code:
12 ENTER 13 RCL+ ST Y

can be replaced by

Code:
12 ENTER 13 STO+ ST X

and

Code:
16 X<>Y 17 Rv 18 + 19 R^

by

Code:
16 STO+ ST Z 17 Rv

saving two bytes! ;-)
Cheers, Werner
09-30-2020, 09:12 AM
Post: #27 Werner Senior Member Posts: 748 Joined: Dec 2013
RE: Lambert W Function (hp-42s)
(09-28-2020 10:01 PM)Albert Chan Wrote:  With Free42, tried x = -1/e:

1 [+/-] [e^X] [+/-] XEQ "eW" ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 0.367879441171

It worked, but result "only" have 17 good digits.
That is a direct consequence of using y'=y'+(y'-y)^2 as a stopping criterion..
You should do 1 extra Newton-Rhapson step
Werner
09-30-2020, 11:04 AM
Post: #28
 lyuka Junior Member Posts: 48 Joined: May 2014
RE: Lambert W Function (hp-42s)
Thanks to Werner, I confirmed that the code was reduced by 2 bytes and updated the web page :-)
09-30-2020, 05:28 PM
Post: #29
 Albert Chan Senior Member Posts: 2,103 Joined: Jul 2018
RE: Lambert W Function (hp-42s)
(09-30-2020 09:12 AM)Werner Wrote:  That is a direct consequence of using y'=y'+(y'-y)^2 as a stopping criterion..
You should do 1 extra Newton-Rhapson step

If y' = y' + (y'-y)^2, it implied y has at least half-precision.
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 half-precision slope. Also, y ln(y) - x = 0. Newton steps does nothing ...
09-30-2020, 05:40 PM (This post was last modified: 09-30-2020 05:42 PM by Albert Chan.)
Post: #30
 Albert Chan Senior Member Posts: 2,103 Joined: Jul 2018
RE: Lambert W Function (hp-42s)
Hi, lyuka

I extended your code to handle complex number input.

Code:
00 { 57-Byte Prgm } 01▶LBL "eW" 02 ENTER 03 ENTER 04 +/- 05 -1 06 E↑X          # r     -x      x      x 07 X=Y? 08 RTN 09 STO- ST Y    # r     -x-r    x      x 10 X<>Y 11 +/-          # x+r   r       x      x 12 ENTER 13 STO+ ST X    # 2x+2r x+r     r      x 14 RCL× ST Z 15 SQRT         # sqrt  x+r     r      x 16 STO+ ST Z 17 R↓           # x+r   sqrt+r  x   sqrt 18 0.3 19 × 20 +            # y     x     x       x 21▶LBL 01 22 STO ST Y 23 LN 24 1 25 + 26 R↑ 27 RCL+ ST Z 28 X<>Y 29 ÷            # Newton-Raphson method 30 - 31 STO× ST X 32 LASTX 33 STO+ ST Y 34 X≠Y?         # Convergence check 35 GTO 01 36 END

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.e-34 - 1.e-33i
09-30-2020, 07:16 PM
Post: #31
 lyuka Junior Member Posts: 48 Joined: May 2014
RE: Lambert W Function (hp-42s)
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 :-)
10-01-2020, 09:37 AM
Post: #32 Werner Senior Member Posts: 748 Joined: Dec 2013
RE: Lambert W Function (hp-42s)
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 @ (Re,0) -> Re @ (Re,Im) -> (Re,Im)  REAL?  0  CPX?  COMPLEX  X#0?  COMPLEX  REAL?  +

or

Code:
 REAL?  GTO 00  COMPLEX  X#0?  COMPLEX  REAL?  +  LBL 00

Cheers, Werner
10-01-2020, 12:33 PM (This post was last modified: 10-01-2020 01:06 PM by Albert Chan.)
Post: #33
 Albert Chan Senior Member Posts: 2,103 Joined: Jul 2018
RE: Lambert W Function (hp-42s)
(10-01-2020 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.158028306006248941790250554076922e-34

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" 02 -1 03 E↑X 04 ENTER 05 RCL+ ST Z    # x+r   r   x 06 PI 07 SIN 08 -            # x+r+ε r   x   x

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.
10-01-2020, 01:39 PM (This post was last modified: 10-01-2020 01:47 PM by Werner.)
Post: #34 Werner Senior Member Posts: 748 Joined: Dec 2013
RE: Lambert W Function (hp-42s)
(10-01-2020 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
10-01-2020, 06:25 PM
Post: #35
 lyuka Junior Member Posts: 48 Joined: May 2014
RE: Lambert W Function (hp-42s)
Using sin(PI) is machine dependent and will take some time.
How about adding +ULP (unit in last place) to the initial value?

Code:
 12 0.3       # 0.3           r+x           r+√(2r(r+x))  x 13 ×         # 0.3(r+x)      r+√(2r(r+x))  x             x 14 LASTX     # 0.3           0.3(r+x)      r+√(2r(r+x))  x 15 1/X       # 1/0.3         0.3(r+x)      r+√(2r(r+x))  x 16 RCL× ST L # 1-ε           0.3(r+x)      r+√(2r(r+x))  x 17 -         # 0.3(r+x)-1+ε  r+√(2r(r+x))  x             x 18 1         # 1             0.3(r+x)-1+ε  r+√(2r(r+x))  x 19 +         # 0.3(r+x)+ε    r+√(2r(r+x))  x             x 20 +         # y0=√(2r(r+x))+0.3(r+x)+ε

Web page updated.
http://www.finetune.co.jp/~lyuka/technot...w-42s.html
10-01-2020, 10:25 PM
Post: #36
 Albert Chan Senior Member Posts: 2,103 Joined: Jul 2018
RE: Lambert W Function (hp-42s)
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.
10-02-2020, 05:44 AM
Post: #37
 lyuka Junior Member Posts: 48 Joined: May 2014
RE: Lambert W Function (hp-42s)
Hi, Albert

Confirmed. And 'eW' on my page is reverted to the previous version 1.1.
I was jumping to the conclusion. Thanks.
10-02-2020, 03:02 PM
Post: #38 Werner Senior Member Posts: 748 Joined: Dec 2013
RE: Lambert W Function (hp-42s)
(09-30-2020 05:28 PM)Albert Chan Wrote:  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.

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*(y-1/e)) + 1 = -1 + LN1P(e*(y-1/e)) + 1
= LN1P(e*(y-1/e))
eg. y = 1/e + 1e-17
then LN(y) + 1 = 2.71828182845904521e-17
LN1P(e*(y-1/e)) = 2.718281828459045198415006976699411e-17

But the cancellation happens also in y + x, and there's nothing we can do. Or at least, nothing *I* can do ;-)

Cheers, Werner
10-02-2020, 05:55 PM (This post was last modified: 10-03-2020 10:24 AM by Albert Chan.)
Post: #39
 Albert Chan Senior Member Posts: 2,103 Joined: Jul 2018
RE: Lambert W Function (hp-42s)
(10-02-2020 03:02 PM)Werner Wrote:  eg. y = 1/e + 1e-17
then ﻿ ﻿ LN(y) + 1 = 2.71828182845904521e-17
LN1P(e*(y-1/e)) = 2.718281828459045198415006976699411e-17

Above amazing accurate results are misleading.
You have y-1/e = 1e-17 (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] 1e-17 [+] ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ // 3.678794411714423315955237701614609e-1 = y

To get an accurate slope, log(y)+1 = log1p(ε = e*y-1) = log1p(ε = e*(y-1/e))
But, this shift the problem to get accurate ε, which required more precise 1/e

1e-17 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ // y - 1/e
-3.255418886896823216549216319830254E-35 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ // (more precise 1/e) - 1/e
− ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ // 1.000000000000000003255418886896823e-17 = y - (more precise 1/e)
1 [EXP] × ﻿ ﻿ ﻿ ﻿ // 2.718281828459045244209433475626668e-17 = ε
[LN1+X] ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ // 2.718281828459045207264152980973417e-17 = log1p(ε)

Mathematica gives 2.718281828459045207264152980973418e-17

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)
10-02-2020, 07:29 PM
Post: #40
 Albert Chan Senior Member Posts: 2,103 Joined: Jul 2018
RE: Lambert W Function (hp-42s)
Here is a demo of previous post results, r + err_r = 1/e (108-bits accurate)

Code:
r, err_r = 0.36787944117144233, -1.2428753672788363E-17 g0 = lambda y,a: (y+a)/(log(y)+1) g1 = lambda y,a: (y+a)/log1p(e*(y - r)) g2 = lambda y,a: (y+a)/log1p(e*(y - r - err_r)) def eW(x, g, n=5):     y = 0.3*(x+r) + sqrt(2*r*(x+r)) + r     if y == r: y += 1e-10   # avoid 0 slope     for i in range(n): print i,y; y = g(y,x)     return y

>>> 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: 1 Guest(s)