[HP41] Lambert function RPN; question

[floppy]

Hello,
I was looking at the forum and internet for finding a lambert function in RPN (my HP41CV dont digest any ClonixD with SANDMATH since yesterday..). I was ready to upload a recommendable function from there in order to understand if/how it works:
https://github.com/isene/hp-41_wlambert/...lambert.41
BUT: What is the line 37 (=RCL M) doing there?
I suppose the 2 programs are calculating L(x) in the high-branch area and low branch area of the lambert calculation. https://en.wikipedia.org/wiki/Lambert_W_function
For the reason the use inputs/outputs are still a bit cryptic for me for the LW and LWM, any remarks/comments are welcome.
Why lambert at Christmas? because "Tetration" is linked to Lambert https://en.wikipedia.org/wiki/Tetration. A RPN code of tetration is included (from Christmas activity..). https://www.desmos.com/calculator/mmhriitkxe?lang=de

Edit 1: link to graphical representation in desmos with the tetration function between X="minimum" and infinite (from e, it is indicated so far with calculation points) was included.

Edit 2: attachment file updated for making it more clear what is the purpose of the function (no programming line updated)

[Didier Lachieze]

(12-26-2022 12:13 PM)floppy Wrote:  I was ready to upload a recommendable function from there in order to understand if/how it works:
https://github.com/isene/hp-41_wlambert/...lambert.41
BUT: What is the line 37 (=RCL M) doing there?
I suppose the 2 programs are calculating L(x) in the high-branch area and low branch area of the lambert calculation. https://en.wikipedia.org/wiki/Lambert_W_function
For the reason the use inputs/outputs are still a bit cryptic for me for the LW and LWM, any remarks/comments are welcome.

Here is the original discussion from where this program is originated : Help solving equation.

You should find there more details about the inputs/outputs of the functions.

Line 37 (=RCL M) is between the end of the LW (RTN) and the beginning of LWM so I suppose it’s there to calculate LWM(x) just after LW(x) by simply pressing R/S.

[Albert Chan]

Two years ago, we had developed a more sophisticated lambertw function (actually e^W)

x*e^x = a      // x = W(a), by definition

Let y = ln(x) --> f = y*ln(y) - a = 0

Newton's method for y: y = y - f/f' = (y+a)/(ln(y)+1)

It extended W for complex numbers, and can get both branches for real argument.

(10-21-2020 03:05 PM)Albert Chan Wrote:  Example, to calculate eW(-log(2)/2), both branch 0 and -1:

2 LN 2 +/- /   → a = -0.3465735902799726547086160607290883
XEQ "eW"       → e^W₀(a) = 0.5
− 0.3 R/S      → e^W_-1(a) = 0.25

For e^W_-1(a), -1/e < a < 0, any positive guess below -a work.

Update: a good guess for e^W_-1(a), -1/e < a < 0, is 2a²
For above case, 2a² already gives 0.24022

Just replace "− new-guess R/S", with this: "− − + × R/S"

^∞c = W(-ln(c)) / -ln(c) = 1 / e^W(-ln(c))

Above quoted example, ^∞(√2) = 2

---

We can reuse infinite tetration formula to solve for a^x = b*x
Let z = b*x, c = ^b√a

z = a^(z/b) = c^z = c^c^z = ... = ^∞c

https://github.com/isene/hp-41_wlambert example, 2^x = 3*x, for x

2 LN 3 / +/-      // -ln(c) = -0.2310490601866484364724107071527255
XEQ "eW"         // 1/z(0) = 0.7280844118213892196256653246326560
− − + × R/S    // 1/z(-1) = 0.1006083268252766116679239025157516

x = z/3

x(0) = 0.4578223732320550555738866680640553
x(-1) = 3.313178380475634845996561019588785

[Sylvain Cote]

(12-26-2022 12:13 PM)floppy Wrote:  my HP41CV dont digest any ClonixD with SANDMATH since yesterday.

Most of Ángel released ROMs are 41CX only. (his ROMs calls 41CX OS sub-routines)

[floppy]

(12-26-2022 02:07 PM)Didier Lachieze Wrote:  Here is the original discussion from where this program is originated : Help solving equation.

You should find there more details about the inputs/outputs of the functions.

Line 37 (=RCL M) is between the end of the LW (RTN) and the beginning of LWM so I suppose it’s there to calculate LWM(x) just after LW(x) by simply pressing R/S.

Thanks. Make sense especially for a debugging.
I moved in-between to a HP41CX and confirm both script for infinite tetration have the same result (that was the expectation):
- above RPN (see attachment) with series; called "TETS" a bit like "TETration infinite Series" (takes.. long..)
- below with lambert WL0 from SANDMATH; called "TETL" a bit like "TETration infinite Lambert"

Quote: 01*LBL "TETL"
02 ENTER^
03 1/X
04 Y^X
05 LN
06 CHS
07 WL0
08 LASTX
09 /
10 RTN
11 .END.

now I will make it 100% RPN according the links above just for my HP41CV
(that diva dont switch on with any clonix/nov module since yesterday; I suppose it dont like any climate change around the house.. but this is another story)
in order to work with standard functions (probably no need to post/update anything more here since the remarks and hints of all made all clear).

THANKS THANKS.

Update: full programm and documentation in TXT attachment

[floppy]

Finally.

Full RPN program for infinite tetration calculation of the parameter value of type "X^(1/X)"; using "LW" 1:1 from the above links; R00 necessary because "LW" screw up the stack (I know, thats not a big program.. but math-prototyping on an HP41 is always fun):

Code:

 01 LBL "TETR"
 02 ENTER^
 03 1/X
 04 Y^X
 05 LN
 06 CHS
 07 STO 00
 08 XEQ "LW"   ; function from this forum and documented on github
 09 RCL 00
 10 /
 11 RTN
 12 END

[Thomas Klemm]

(12-26-2022 02:07 PM)Didier Lachieze Wrote:  Here is the original discussion from where this program is originated : Help solving equation.

It took me a while to figure out that I had used Halley's method with:

\(
\begin{align}
f(x) &= x e^x - a \\
f'(x) &= (1 + x) e^x \\
f''(x) &= (2 + x) e^x \\
\end{align}
\)

This is the formula used in the program:

\(
\begin{align}
x_{{n+1}}
&=x_{n}-{\frac {f(x_{n})}{f'(x_{n})-{\frac {f(x_{n})}{f'(x_{n})}}{\frac {f''(x_{n})}{2}}}} \\
\\
&=x_{n}-{\frac{x e^x - a}{(1 + x) e^x-\frac{(x e^x - a)(2 + x) e^x}{2(1 + x) e^x}}} \\
\\
&=x_{n}-{\frac{x e^x - a}{(1 + x) e^x-\frac{(x e^x - a)(2 + x)}{2(1 + x)}}} \\
\\
&=x_{n}-{\frac{x e^x - a}{(1 + x) e^x-\frac{(x e^x - a)(1 + \frac{1}{1 + x})}{2}}} \\
\end{align}
\)

But we could also use:

\(
\begin{align}
x_{n+1}
&=x_{n}-{\frac {f(x_{n})}{f'(x_{n})}}\left[1-{\frac {f(x_{n})}{f'(x_{n})}}\cdot \frac {f''(x_{n})}{2f'(x_{n})}\right]^{-1} \\
\\
&=x_{n}-{\frac {x e^x - a}{(1 + x) e^x}}\left[1-{\frac {x e^x - a}{(1 + x) e^x}}\cdot \frac {(2 + x) e^x}{2(1 + x) e^x}\right]^{-1} \\
\\
&=x_{n}+{\frac {a - x e^x}{e^x + x e^x}}\left[1+{\frac {a - x e^x}{e^x + x e^x}}\cdot \frac {2 + x}{2 + 2x}\right]^{-1} \\
\end{align}
\)

We end up with a slightly shorter program, here for the HP-42S:

Code:

00 { 41-Byte Prgm }
01▸LBL "W"
02 STO 00
03 0
04 STO 01
05▸LBL 00
06 RCL 00
07 RCL 01
08 ENTER
09 E↑X
10 ×
11 STO- ST Y
12 LASTX
13 +
14 ÷
15 RCL 01
16 2
17 RCL 01
18 +
19 +
20 LASTX
21 ÷
22 ÷
23 1
24 +
25 ÷
26 STO+ 01
27 ABS
28 RND
29 X≠0?
30 GTO 00
31 RCL 01
32 END

But it should work for the HP-41C as well.

Example

4 XEQ "W"

1.20216787320

So why Halley and not Newton, you might ask.
With a little more effort, we get cubic convergence, which is nice.

[Albert Chan]

(12-27-2022 02:31 AM)Thomas Klemm Wrote:  It took me a while to figure out that I had used Halley's method with:

\(
\begin{align}
f(x) &= x e^x - a \\
f'(x) &= (1 + x) e^x \\
f''(x) &= (2 + x) e^x \\
\end{align}
\)

There was a problem with this setup. Derivatives may be huge, possibly bigger than f itself!
Had we use simple Newton's method, it may not converge at all.

>>> from mpmath import *
>>> a = 1e99
>>> x = ln(a) # W(a) guess ≈ 228
>>> def halley(x): y=exp(x); z=x*y; f=z-a; return f/(y+z - f/(y+z)*(y+z/2))
...
>>> for i in range(6): x -= halley(x); print x
...
225.982091142775
224.113871381155
222.808520559707
222.552199206554
222.550768955996
222.55076895575

Solve x*exp(x) = a, Halley's method, take 6 iterations.
Solve y*ln(y) = a, y=e^x, Newton's method only take 4.

>>> y = a # e^W(a) guess
>>> for i in range(4): y=(y+a)/(log(y)+1); print a/y
...
114.477962103205
222.273911298397
222.550768184143
222.55076895575

A better Halley setup is first temper down huge derivatives.
Bonus: Halley's correction term is simpler.

f = x + ln(x/a)
f' = 1 + 1/x
f'' = -1/x^2

>>> x = ln(a) # same guess
>>> def halley(x): f=x+ln(x/a); return f*x/(x+1 + f/(x+x+2))
...
>>> for i in range(2): x -= halley(x); print x
...
222.550764466154
222.55076895575

[Thomas Klemm]

Nice. This removes another 4 bytes:

Code:

00 { 37-Byte Prgm }
01▸LBL "W"
02 STO 00
03 LN
04 STO 01
05▸LBL 00
06 RCL 01
07 RCL 01
08 RCL 00
09 ÷
10 LN
11 RCL 01
12 +
13 ×
14 LASTX
15 2
16 ÷
17 RCL 01
18 1
19 +
20 ÷
21 LASTX
22 +
23 ÷
24 STO- 01
25 ABS
26 RND
27 X≠0?
28 GTO 00
29 RCL 01
30 END

In addition, it is significantly faster.

Example

1E99 XEQ "W"

222.550768956

[Thomas Klemm]

(12-26-2022 06:45 PM)floppy Wrote:  R00 necessary because "LW" screw up the stack

Since "LW" saves the argument in register M you could also just recall it from there:

Code:

 01 LBL "TETR"
 02 ENTER^
 03 1/X
 04 Y^X
 05 LN
 06 CHS
 07 XEQ "LW"   ; function from this forum and documented on github
 08 RCL M
 09 /
 10 END

The END command serves also as RTN.

[Thomas Klemm]

(12-27-2022 04:26 AM)Albert Chan Wrote:  print x

Are you seriously still using Python 2?
This parrot is dead.
Sunsetting Python 2

[floppy]

(12-27-2022 04:48 PM)Thomas Klemm Wrote:  

Code:

 01 LBL "TETR"
 02 ENTER^
 03 1/X
 04 Y^X
 05 LN
 06 CHS
 07 XEQ "LW"   ; function from this forum and documented on github Geir Isene
 08 RCL M
 09 /
 10 END

+1

[Pekis]

Here's a simple algorithm derived from the one I found on math.stackexchange for both branches of the Lambert function:

The iteration
Using Newton's method to solve x * exp(x) = p yields the following iterative step for finding x = W(p):
x' = (p * exp(−x) + x * x) / (x + 1)

Initial values x0
For the W0 branch,
When p is in [−1/e, 10], use x0 = 0
When p>10, use x0 = ln(p)−ln(ln(p)), ie following asymptotic behavior

For the W-1 branch,
When p is in [−1/e,−0.1], use x0 = −2
When p is in ]−0.1,0], use x0 = ln(−p)−ln(−ln(−p)), ie following asymptotic behavior

[Albert Chan]

(04-10-2022 04:03 AM)Albert Chan Wrote:  I borrowed lyuka e^W estimation formula, replace 0.3 by 1/3, for Puiseux series.

e^W(y) ≈ 1/e + sqrt ((2/e)*(y+1/e)) + (y+1/e)/3

If we replace +sqrt (for 0 branch), with -sqrt (for -1 branch), formula also work (good enough for guess)

Formula is exact where 2 branches intersect, e^W(-1/e) = 1/e = 0.367879 441171 ...

e^W₀ (-0.367879) ≈ 0.368449 321330
e^W_-1(-0.367879) ≈ 0.367309 855127
e^W₀ (-0.3678) ≈ 0.375551 151939
e^W_-1(-0.3678) ≈ 0.360260 691185

Update: unfortunately, it does not work well for argument close to 0
Branch -1 with argument clos to 0, e^W guess may return negative number.
Too bad ...

[Albert Chan]

When I try to understand lyuka's eW formula (previous post), I found a better one.

e^W(a) = y = (y+a) / (ln(y)+1)

Let r = 1/e
If a = -r + h (tiny h), we have y = r + z (tiny z)

e^W(a) = (r+z) = (h+z) / ln(1+e*z)

Let H = e*h (known), Z = e*z (unknown), we have:

(1+Z) = (H+Z) / ln(1+Z)

H = (1+Z) * ln(1+Z) - Z = Z^2/2 - Z^3/6 + Z^4/12 - Z^5/20 + ...

2H ≈ Z^2 / (1+Z/3)
Z^2 - 2*(H/3)*Z - 2H ≈ 0
Z ≈ H/3 ± √((H/3)^2 + 2H) ≈ H/3 ± √(2H)

Divide both side by e, we have: y = r + z ≈ r + h/3 ± √(2*r*h)

We could estimate Z another way, by first ignoring O(Z^3) terms
Note: +sign for branch 0, -sign for branch -1

2H ≈ Z^2                      → Z = ±√(2H)

Then, we solve again, using above rough estimated Z

2H ≈ Z^2 / (1+Z/3)      → Z = ±√(2H*(1+Z/3))

Redo previous examples, with this 2 step Z estimate, then convert back to y = r + r*Z

e^W₀ (-0.367879) ≈ 0.368449 321311
e^W_-1(-0.367879) ≈ 0.367309 855146
e^W₀ (-0.3678) ≈ 0.375551 106237
e^W_-1(-0.3678) ≈ 0.360260 737204

[Albert Chan]

The reason for accurate eW estimate at a ≈ -1/e, is that it take a lot of iterations for convergence.

Once we get close, convergence speed up dramatically (slope = ln(1+Z) ≈ Z)
For IEEE double, if we get around 40 bits, Newton's step converged in 1 step.

a ≈ -1/e --> y ≈ 1/e, and both are inputs (considered exact) --> (y+a) is exact.

We do need to use log1p() instead of ln() for accurate slope.

ln(y)+1 = ln(e*y) = log1p(e*(y-1/e))

Term (1/e) need extra precisions, to avoid cancellation errors.
--> (r + err) = rounded 108 bits of 1/e (about 33 decimals digits accuracy)

Code:

function expW(a, y, verbal)
    y = y or 0  -- default branch 0
    local h, s
    if a <= -0.1 then
        h = a + r + err
        y = sqrt(2*r*h) * (y+y+1)
        y = y * sqrt(1+y/(3*r)) + r
        s = function(y) return log1p((y-r-err)/r) end
    else
        y = (y+1) + (y+y+1)*a/4
        s = function(y) return log(y) + 1 end
    end
    repeat
        if verbal then print(y, a/y) end
        h = y - (y+a) / s(y)
        y = y - h   -- Newton's correction
    until y == y+h*0x1p-13 or y ~= y
    return y, a/y   -- e^W(a), W(a)
end

Convergence should take at most 6 iterations
Accuracy should be +/- a few ulps, for a ≥ -1/e

lua> expW(1e300) -- default = branch 0
1.4614601088436296e+297      684.2472086297608
lua> _ * log(_)
1e+300

lua> expW(-1e-300, -1) -- branch -1
1.4340561272249246e-303      -697.3227762954601
lua> _ * log(_)
-1e-300

Update 1: I implemented Lua code for x = W(a), solving for f = x + ln(x/a) = 0
This post is updated to match that, with a verbal option.
Also, guard removed, "if h == err then return r, -1 end"

Update 2: change behavior from "W branch -1 if y" to "y default to branch 0"
Assumed input, y = nil/false/0/-1 only: y = nil/false/0 get branch 0; y = -1 get branch -1

lua> expW(-0.367, 0, true)
0.3936082265591097              -0.9323992112875368
0.3936082377626915              -0.9323991847479225
0.3936082377626891              -0.9323991847479282