HP Forums - Lambert W Function (hp-42s)

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

What's wrong with this?

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

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?

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

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

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 Smile

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