# HP Forums

Full Version: Lambert W Function (hp-42s)
You're currently viewing a stripped down version of our content. View the full version with proper formatting.
Pages: 1 2 3
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.
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
(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
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
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
(05-17-2020 09:29 AM)Gerald H Wrote: [ -> ]What's wrong with this?

It has two global labels ;-)
Werner
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
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.
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
(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
Hello Albert, how very interesting your equation, and it looks simpler too.
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
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
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
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
Hi Albert,
something must be missing? Your routine expects two arguments.
(the '+' in line 3)
Werner
Hi, Werner

Thanks. Code corrected.

BTW, can the code be done all in the stack (no "a") ?
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```
(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
Pages: 1 2 3
Reference URL's
• HP Forums: https://www.hpmuseum.org/forum/index.php
• :