Lambert W Function (hp-42s)

[Juan14]

The Lambert W Function can be implemented in your calculator using Newton's Method. To find W(A), the equation to solve is:
f(x) = x*exp(x)-A = 0 (1)
f”(x) = (x+1)*exp(x) (2)

A recursive formula for Newton's Method is:
X = x+(A/exp(x)-x)/(x+1)) (3)

From (1) we have x*exp(x) = A or ln(x) + x = ln(A); x = ln(A)-ln(x); x ~ ln(A); ln(A) is a good initial value for the Newton Method when A>1, ln(A+1.4) is a good initial value for all reals ≥ -1/e, for complex numbers we use ln(A).

The next program implements the Lambert W Function For the hp 42s:

01 LBL “LWF”
02 1.4
03 X<>Y
04 ENTER
05 REAL?
06 RCL+ ST Z
07 LN
08 STO “X”
09 LBL 00
10 X<>Y
11 RCL ST X
12 RCL “X”
13 E↑X
14 ÷
15 RCL- “X”
16 1
17 RCL+ “X”
18 ÷
19 STO+ “X”
20 ABS
21 1E-32
22 -
23 X>0?
24 GTO 00
25 RCL “X”
26 END

The program doesn't check if you entered a value where the function its not defined. Values of x < -1/e should be entered as complex.
The strategy to solve equations with exponential using the Lambert W Function is to use algebraic manipulations to get an equation of the form x*exp(x) = a; then x = W(a)
or x*ln(x) = a; then ln(x) = W(a) and x = exp(W(a)).

Example: to solve x^x = π, we have ln(x^x) = ln(π), x*ln(x) = ln(π); x = exp(W(ln(π))); in the calculator:

[■][π][LN][XEQ][LWF][■][e^x] give us a value of x ~ 1.85410596792.

[Werner]

Why not use the built-in solver?

Code:

00 { 36-Byte Prgm }
01*LBL "W"
02 FC? 45
03 GTO 00
04 RCL "X"
05 E^X
06 LASTX
07 x
08 RCL- "W"
09 RTN
10*LBL 00
11 PGMSLV "W"
12 LSTO "W"
13 CLX
14 LSTO "X"
15 10
16 SOLVE "X"
17 END

Cheers, Werner

[J-F Garnier]

(05-17-2020 07:56 AM)Werner Wrote:  

Code:

00 { 36-Byte Prgm }
01*LBL "W"
02 FC? 45
03 GTO 00
...
10*LBL 00
11 PGMSLV "W"
12 LSTO "W"
13 CLX
14 LSTO "X"
15 10
16 SOLVE "X"
17 END

Nice use of flag 45.

Of course, LSTO doesn't exist on the HP-42S.

J-F

[Werner]

I'm afraid Free42/DM42 has taken over for me ;-)
(actually, the version that I have stored is with STO; but I wanted to see if it worked with LSTO, and of course it does. I will not doubt you again, Thomas.)

Cheers, Werner

[Gerald H]

What's wrong with this?

https://www.hpmuseum.org/forum/thread-76...ht=Lambert

[Juan14]

I don't own a hp-42s, line 21 should be 1E-10 or something like that in the hp-42s. My program can take advantage of the fact that the hp-42s can handle complex numbers. For example let's consider a infinite tower of the imaginary number i, i^i^i^i......
We have i^i^i^i...... = x; i^(i^i^i^i)...... = x; i^x = x or i = x^(1/x) taking ln on both sides of the equation
ln(i) = 1/x ln(x); -ln(i) = -1/x ln(x); -ln(i) = 1/x ln(1/x), calling W to the Lambert W function.
W(-ln(i)) = W(1/x ln(1/x)); W(-ln(i)) = ln(1/x); 1/x = exp(W(-ln(i))) and finally x=exp(-W(-ln(i))).
In the calculator:
0 [ENTER] 1 [COMPLEX] [LN] [+/-] [XEQ] LWF [+/-][■][E↑X]
You get 4.3828293673E-1 i 3.605924718714E-1
so i^i^i^i...... ≈ 0.43828293673+ i 0.3605924718714

[Werner]

(05-17-2020 09:29 AM)Gerald H Wrote:  What's wrong with this?

https://www.hpmuseum.org/forum/thread-76...ht=Lambert

It has two global labels ;-)
Werner

[Albert Chan]

A better convergence test might be (x-prev_x)*eps + x - x = 0

With Newton's method, eps can be set to "half" the calculator precision.
Example, 12 digits calculator, eps = 1e-6 is sufficient.

Quote:For example let's consider a infinite tower of the imaginary number i, i^i^i^i......

Instead of using LambertW, we may also do the direct route.

x = k^k^k ... = k^x
log(x) = x log(k)
Let f = log(x) - x log(k), f' = 1/x - log(k)

Newton's iteration x = x - f(x)/f'(x):

x = (1 - log(x)) / (1/x - log(k)
y = 1/x = (y - log(k)) / (log(y) + 1)

>>> from cmath import *
>>> a = -log(1j)                                     # k = 1j
>>> g = lambda y,a: (y+a)/(log(y)+1)     # g(1,a) = 1+a
>>> g(1+a,a) →(1.3127784610672455-1.1245675977983851j)
>>> g(_,a)     →(1.3607126098119278-1.119058042805926j)
>>> g(_,a)     →(1.3606248500898628-1.1194391815927829j)
>>> g(_,a)     →(1.3606248702911174-1.1194391662423497j)
>>> g(_,a)     →(1.3606248702911177-1.11943916624235j), pass convergence test, eps=2**-26
>>> y = _
>>> x = 1/y   →(0.43828293672703211+0.36059247187138549j)
>>> 1j ** x    →(0.43828293672703211+0.36059247187138554j)

x = 1/y = W(a)/a → (a/y) exp(a/y) = a → exp(a/y) = y → W(a) = a/y = log(y)

>>> a/y         → (0.56641733028546448-0.68845322710770207j)
>>> log(y)     → (0.56641733028546437-0.68845322710770218j)

Comment: Solving for x = k^x is not quite the same as x = k^k^k ...
Example:

3 = (³√3)^3, but (³√3) ^ (³√3) ^ (³√3) ... ≈ 2.47805268029 ≠ 3

[Juan14]

Hello Albert, I know that any equation that can be solved by Lambert W Function, can be solved using a numeric method, but I found satisfying using the algebraic manipulation so I can use Lambert W Function.

[Albert Chan]

Hi, Juan14

I was experimenting another way to get W(a), via infinite tetration route.
Solving for y = exp(W(a)) maybe better than directly going for W(a):

Code:

from cmath import *
g = lambda y,a: (y+a)/(log(y)+1)
w = lambda x,a: x+(a/exp(x)-x)/(x+1)

def fixedpoint(f, x, a, eps=2**-26, limit=100):
    for t in xrange(1,limit):   # f() calls
        prev, x = x, f(x,a)
        if (x-prev)*eps + x == x: return x, t

def test(a):
    y, t = fixedpoint(g, 1+a, a)
    print 'y->w =', (a/y, t)
    print 'w    =', fixedpoint(w, log(1+a), a)

>>> test(1)
y->w = ((0.56714329040978384+0j), 4)
w      = ((0.56714329040978384+0j), 5)
>>> test(1e10)
y->w = ((20.028685413304952+0j), 5)
w      = ((20.028685413304952+0j), 8)
>>> test(1e100)
y->w = ((224.84310644511851+0j), 4)
w      = ((224.84310644511851+0j), 10)

Note: both versions use the same guess, y = e^x = 1+a
Note: y->w uses a/y instead of log(y) to recover W(a), to reduce errors when y ≈ 1

[Albert Chan]

(05-19-2020 12:29 AM)Albert Chan Wrote:  g = lambda y,a: (y+a)/(log(y)+1)
w = lambda x,a: x+(a/exp(x)-x)/(x+1)

Both Newton's iteration formula are equivalent, but going for y converge faster.
(try plotting x*exp(x) vs x*ln(x), 2nd curve slope = 1+ln(x), does not vary as much)

W(a) = x
x * exp(x) = a

Let y = exp(x), we have y * ln(y) = a

Newton's iteration: y ← y - (y*ln(y)-a) / (ln(y)+1) = (y+a) / (ln(y)+1)

When y converged, y = exp(x) = a/x, we have W(a) = ln(y) = a/y

[Juan14]

Hello Albert, how very interesting your equation, and it looks simpler too.

[Werner]

Ah yes, my simple program of course only works for reals.
Okay then, let's use a general purpose complex solver.
I converted the one written by Valentin Albillo for the HP-35S, to Free42/DM42.
(you can find it here on his site)
Of course, it is equally usable on the 42S if you change all LSTO statements to STO, and change the 1E-15 to 1E-5.
I changed three things in the program:
- a more accurate version of the root of the quadratic
- a simpler test to see whether the number is sufficiently close to a real (that is not possible on a 35S or Valentin would've used it)
- and just taking the real part of the current estimate if the test is positive

Usage:
ALPHA: name of function to solve
X: guess
XEQ "CSLV"

The guess may be real or complex, and the resulting root as well.

Code:

00 { 145-Byte Prgm }
01▸LBL 01
02 RCL+ "X"
03 ASTO ST L
04 GTO IND ST L
05▸LBL "CSLV"
06 0
07 LSTO "U"
08 LSTO "V"
09 ENTER
10 COMPLEX
11 +
12 LSTO "X"
13 1ᴇ-15
14 LSTO "S"
15▸LBL 02
16 CLX
17 XEQ 01
18 STO+ ST X
19 STO "U"
20 RCL "S"
21 +/-
22 XEQ 01
23 STO "V"
24 RCL "S"
25 XEQ 01
26 ENTER
27 RCL+ "V"
28 RCL- "U"
29 RCL "S"
30 X↑2
31 ÷
32 X<>Y
33 RCL- "V"
34 RCL "S"
35 STO+ ST X
36 ÷
37 ÷
38 RCL "U"
39 LASTX
40 ÷
41 STO ST Z
42 ×
43 1
44 -
45 +/-
46 SQRT
47 1
48 +
49 ÷
50 STO- "X"
51 RCL÷ "X"
52 ABS
53 RCL "S"
54 X↑2
55 X<Y?
56 GTO 02
57 RCL "X"
58 SIGN
59 COMPLEX
60 X<>Y
61 R↓
62 ABS
63 X>Y?
64 GTO 00
65 RCL "X"
66 COMPLEX
67 R↓
68 STO "X"
69▸LBL 00
70 RCL "X"
71 END

So, to find x=i^i^i..., solve f(x) = i^x - x = 0

define FX:

Code:

>LBL "FX"
 -1
 SQRT
 X<>Y
 Y^X
 LASTX
 -
 END

"FX"
0
XEQ "CSLV"
(0.43829..,0.36059..)
XEQ "FX" results in (0,-1e-34)

Implement Lambert's W:

Code:

>LBL "W"
 LSTO "W"
 "WX"
 1.4
 X<>Y
 REAL?
 +
 XEQ "CSLV"
 RTN
>LBL "WX"
 E^X
 LASTX
 x
 RCL- "W"
 END

The solution to the above problem can be written as

x = i^x

x = e^(x*ln(i))

x*ln(i) = ln(i)*e^(x*ln(i))

ln(i) = x*ln(i)*e^(-x*ln(i))

Let y = -x*ln(i), then solve

y^*e^y = -ln(i) or calculate W(-ln(i))

-1
SQRT
LN
+/-
XEQ "W"
-1
SQRT
LN
+/-
/

results in the almost the the same number as above.

However, Houston, we have a problem here. The definition of WX must not use the local variables used by CSLV, namely X, S, U and V. If we had used variable X instead of W, the routine would not have worked. To remove this dependency, here's a version of CSLV that does not use alphanumeric variables, except REGS:
(making it less usable on a real 42S because it will overwrite REGS)

00 { 115-Byte Prgm }
01▸LBL 01
02 RCL 00
03 +
04 ASTO ST L
05 GTO IND ST L
06▸LBL "CSLV"
07 2
08 ENTER
09 NEWMAT
10 ENTER
11 COMPLEX
12 +
13 LSTO "REGS"
14 1ᴇ-15
15 STO 01
16▸LBL 02
17 CLX
18 XEQ 01
19 STO+ ST X
20 STO 02
21 RCL 01
22 +/-
23 XEQ 01
24 STO 03
25 RCL 01
26 XEQ 01
27 ENTER
28 RCL+ 03
29 RCL 02
30 -
31 RCL 01
32 X↑2
33 ÷
34 X<>Y
35 RCL 03
36 -
37 RCL 01
38 STO+ ST X
39 ÷
40 ÷
41 RCL 02
42 LASTX
43 ÷
44 STO ST Z
45 ×
46 1
47 -
48 +/-
49 SQRT
50 1
51 +
52 ÷
53 STO- 00
54 RCL 00
55 ÷
56 ABS
57 RCL 01
58 ABS
59 X↑2
60 X<Y?
61 GTO 02
62 RCL 00
63 SIGN
64 COMPLEX
65 X<>Y
66 R↓
67 ABS
68 X>Y?
69 GTO 00
70 RCL 00
71 COMPLEX
72 R↓
73 STO 00
74▸LBL 00
75 RCL 00
76 END

Hope you like it, Werner

[Albert Chan]

Hi, Werner

Very nice complex solver.
Regarding LambertW code, it missed LN before entering solver.
Instead of guess of LN(1.4+X) or LN(X), I found that guess LN(1+X) is simpler.

Here is the revised LambertW (34-byte)

Code:

01▸LBL "W"
02 LSTO "W"
03 "WX"
04 1
05 +
06 LN
07 XEQ "CSLV"
08 RTN
09▸LBL "WX"
10 E↑X
11 LASTX
12 ×
13 RCL- "W"
14 END

Example, i^i^i ... = W(-ln(i))/(-ln(i)) = exp(-W(-ln(i)))

-1
SQRT
LN
+/-
XEQ "W"
+/-
EXP

→ 0.4382829367270321116269751635512648 + 0.3605924718713854859529405269060007i

[Werner]

Thanks Albert. You can then simply use LN1+X ;-)
Update: no you can't! LN accepts complex arguments, but LN1+X does not. Another example of the incredible attention to detail in Free42/DM42.
Cheers, Werner

[Albert Chan]

This gives exp(W(a)), which equals a / W(a)

Code:

00 { 33-Byte Prgm }
01▸LBL "eW"
02 LSTO "a"
03 1
04 +
05▸LBL 01
06 ENTER        ; y     y
07 RCL+ "a"     
08 LASTX        ; y     a+y     y
09 LN
10 1
11 +
12 ÷            ; y'    y
13 -            ; eps
14 LASTX        
15 X<>Y         ; eps   y'
16 X↑2
17 X<>Y         ; y'    eps^2
18 +
19 LASTX        ; y'    y'+eps^2
20 X≠Y?
21 GTO 01
22 END

Example, for i^i^i ...

-1
SQRT
LN
+/-
XEQ "eW"
1/X

→ 0.438282936727 + 0.360592471871 j

[Werner]

Hi Albert,
something must be missing? Your routine expects two arguments.
(the '+' in line 3)
Werner

[Albert Chan]

Hi, Werner

Thanks. Code corrected.

BTW, can the code be done all in the stack (no "a") ?

[Werner]

Sure. and if you replace lines 10 and 11 by
X<>Y
STO+ ST Z
X<> ST Z
it is even 41-compatible.

Code:

00 { 33-Byte Prgm }
01▸LBL "eW"
02 STO ST Z
03 1
04 +
05▸LBL 01 @       y               a
06 STO ST Y
07 LN
08 1
09 + @            1+LN(y) y       a       a
10 R↑
11 RCL+ ST Z @    a+y     1+LN(y) y       a
12 X<>Y
13 ÷ @            y'      y       a       a
14 -
15 STO× ST X
16 LASTX
17 STO+ ST Y @    y'    y'+e^2    a    a
18 X≠Y?
19 GTO 01
20 END

[Albert Chan]

(06-11-2020 09:20 AM)Werner Wrote:  

Code:

17 STO+ ST Y @    y'    y'+e^2    a    a

Thanks ! This is more than I expected
Now, we can get W both ways, directly from the stack

1e100
XEQ "eW"   ; eW = 4.44754573894e97
LN         ; W = LN(eW) = 224.843106445
R↓         ; a = 1e100
÷          ; W = a / eW = 224.843106445

Or, to check its accuracy, from definition W eW - a = 0

1e100 XEQ "eW" LN × −     → 0