+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html)
+--- Forum: General Forum (/forum-4.html)
+--- Thread: Lambert W Function (hp-42s) (/thread-15004.html)
Lambert W Function (hp-42s) - Juan14 - 05-16-2020 04:07 PM
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.
RE: Lambert W Function (hp-42s) - Werner - 05-17-2020 07:56 AM
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 ENDCheers, Werner
RE: Lambert W Function (hp-42s) - J-F Garnier - 05-17-2020 08:09 AM
(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
RE: Lambert W Function (hp-42s) - Werner - 05-17-2020 08:15 AM
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
RE: Lambert W Function (hp-42s) - Gerald H - 05-17-2020 09:29 AM
What's wrong with this?
https://www.hpmuseum.org/forum/thread-7692.html?highlight=Lambert
RE: Lambert W Function (hp-42s) - Juan14 - 05-17-2020 12:12 PM
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
RE: Lambert W Function (hp-42s) - Werner - 05-18-2020 08:04 AM
(05-17-2020 09:29 AM)Gerald H Wrote: What's wrong with this?
https://www.hpmuseum.org/forum/thread-7692.html?highlight=Lambert
It has two global labels ;-)
Werner
RE: Lambert W Function (hp-42s) - Albert Chan - 05-18-2020 03:54 PM
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
RE: Lambert W Function (hp-42s) - Juan14 - 05-18-2020 10:51 PM
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.
RE: Lambert W Function (hp-42s) - Albert Chan - 05-19-2020 12:29 AM
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
RE: Lambert W Function (hp-42s) - Albert Chan - 05-20-2020 03:25 PM
(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
RE: Lambert W Function (hp-42s) - Juan14 - 05-21-2020 12:09 AM
Hello Albert, how very interesting your equation, and it looks simpler too.
RE: Lambert W Function (hp-42s) - Werner - 05-22-2020 11:39 AM
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 ENDdefine 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"
ENDx = 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
RE: Lambert W Function (hp-42s) - Albert Chan - 05-23-2020 12:43 AM
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 ENDExample, i^i^i ... = W(-ln(i))/(-ln(i)) = exp(-W(-ln(i)))
-1
SQRT
LN
+/-
XEQ "W"
+/-
EXP
→ 0.4382829367270321116269751635512648 + 0.3605924718713854859529405269060007i
RE: Lambert W Function (hp-42s) - Werner - 05-23-2020 04:20 AM
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
RE: Lambert W Function (hp-42s) - Albert Chan - 06-11-2020 03:27 AM
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 ENDExample, for i^i^i ...
-1
SQRT
LN
+/-
XEQ "eW"
1/X
→ 0.438282936727 + 0.360592471871 j
RE: Lambert W Function (hp-42s) - Werner - 06-11-2020 05:17 AM
Hi Albert,
something must be missing? Your routine expects two arguments.
(the '+' in line 3)
Werner
RE: Lambert W Function (hp-42s) - Albert Chan - 06-11-2020 08:12 AM
Hi, Werner
Thanks. Code corrected.
BTW, can the code be done all in the stack (no "a") ?
RE: Lambert W Function (hp-42s) - Werner - 06-11-2020 09:20 AM
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 ENDRE: Lambert W Function (hp-42s) - Albert Chan - 06-11-2020 10:42 AM
(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