Lambert W Function (hp-42s)
05-16-2020, 04:07 PM
Post: #1
 Juan14 Junior Member Posts: 43 Joined: Jan 2014
Lambert W Function (hp-42s)
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.
05-17-2020, 07:56 AM
Post: #2
 Werner Senior Member Posts: 747 Joined: Dec 2013
RE: Lambert W Function (hp-42s)
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, 08:09 AM
Post: #3
 J-F Garnier Senior Member Posts: 752 Joined: Dec 2013
RE: Lambert W Function (hp-42s)
(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
05-17-2020, 08:15 AM
Post: #4
 Werner Senior Member Posts: 747 Joined: Dec 2013
RE: Lambert W Function (hp-42s)
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
05-17-2020, 09:29 AM
Post: #5
 Gerald H Senior Member Posts: 1,573 Joined: May 2014
RE: Lambert W Function (hp-42s)
What's wrong with this?

05-17-2020, 12:12 PM
Post: #6
 Juan14 Junior Member Posts: 43 Joined: Jan 2014
RE: Lambert W Function (hp-42s)
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-18-2020, 08:04 AM
Post: #7
 Werner Senior Member Posts: 747 Joined: Dec 2013
RE: Lambert W Function (hp-42s)
(05-17-2020 09:29 AM)Gerald H Wrote:  What's wrong with this?

It has two global labels ;-)
Werner
05-18-2020, 03:54 PM (This post was last modified: 05-22-2020 03:08 PM by Albert Chan.)
Post: #8
 Albert Chan Senior Member Posts: 2,103 Joined: Jul 2018
RE: Lambert W Function (hp-42s)
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
05-18-2020, 10:51 PM
Post: #9
 Juan14 Junior Member Posts: 43 Joined: Jan 2014
RE: Lambert W Function (hp-42s)
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.
05-19-2020, 12:29 AM (This post was last modified: 05-25-2020 12:14 PM by Albert Chan.)
Post: #10
 Albert Chan Senior Member Posts: 2,103 Joined: Jul 2018
RE: Lambert W Function (hp-42s)
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-20-2020, 03:25 PM (This post was last modified: 05-20-2020 05:37 PM by Albert Chan.)
Post: #11
 Albert Chan Senior Member Posts: 2,103 Joined: Jul 2018
RE: Lambert W Function (hp-42s)
(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
05-21-2020, 12:09 AM
Post: #12
 Juan14 Junior Member Posts: 43 Joined: Jan 2014
RE: Lambert W Function (hp-42s)
Hello Albert, how very interesting your equation, and it looks simpler too.
05-22-2020, 11:39 AM (This post was last modified: 05-22-2020 11:41 AM by Werner.)
Post: #13
 Werner Senior Member Posts: 747 Joined: Dec 2013
RE: Lambert W Function (hp-42s)
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
05-23-2020, 12:43 AM (This post was last modified: 05-23-2020 12:48 AM by Albert Chan.)
Post: #14
 Albert Chan Senior Member Posts: 2,103 Joined: Jul 2018
RE: Lambert W Function (hp-42s)
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
05-23-2020, 04:20 AM (This post was last modified: 05-23-2020 10:02 AM by Werner.)
Post: #15
 Werner Senior Member Posts: 747 Joined: Dec 2013
RE: Lambert W Function (hp-42s)
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
06-11-2020, 03:27 AM (This post was last modified: 06-11-2020 07:58 AM by Albert Chan.)
Post: #16
 Albert Chan Senior Member Posts: 2,103 Joined: Jul 2018
RE: Lambert W Function (hp-42s)
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
06-11-2020, 05:17 AM
Post: #17
 Werner Senior Member Posts: 747 Joined: Dec 2013
RE: Lambert W Function (hp-42s)
Hi Albert,
something must be missing? Your routine expects two arguments.
(the '+' in line 3)
Werner
06-11-2020, 08:12 AM (This post was last modified: 06-11-2020 08:14 AM by Albert Chan.)
Post: #18
 Albert Chan Senior Member Posts: 2,103 Joined: Jul 2018
RE: Lambert W Function (hp-42s)
Hi, Werner

Thanks. Code corrected.

BTW, can the code be done all in the stack (no "a") ?
06-11-2020, 09:20 AM
Post: #19
 Werner Senior Member Posts: 747 Joined: Dec 2013
RE: Lambert W Function (hp-42s)
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, 10:42 AM
Post: #20
 Albert Chan Senior Member Posts: 2,103 Joined: Jul 2018
RE: Lambert W Function (hp-42s)
(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
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)