Lambert W Function (hp-42s)

[lyuka]

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...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.

[Albert Chan]

(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)

[lyuka]

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.

[Albert Chan]

(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 ...

[lyuka]

That case of x = -1 / e has already been handled by my code shown in the link in the previous post.

[Werner]

(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

[Werner]

(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

[lyuka]

Thanks to Werner, I confirmed that the code was reduced by 2 bytes and updated the web page :-)

[Albert Chan]

(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 ...

[Albert Chan]

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

[lyuka]

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 :-)

[Werner]

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

[Albert Chan]

(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.

[Werner]

(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

[lyuka]

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

[Albert Chan]

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.

[lyuka]

Hi, Albert

Confirmed. And 'eW' on my page is reverted to the previous version 1.1.
I was jumping to the conclusion. Thanks.

[Werner]

(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

[Albert Chan]

(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)

[Albert Chan]

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')