[HP41] Lambert function RPN; question - Printable Version +- 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: [HP41] Lambert function RPN; question (/thread-19342.html) [HP41] Lambert function RPN; question - floppy - 12-26-2022 12:13 PM 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/blob/master/hp41wlambert.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) RE: [HP41] Lambert function RPN; question - Didier Lachieze - 12-26-2022 02:07 PM (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/blob/master/hp41wlambert.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. RE: [HP41] Lambert function RPN; question - Albert Chan - 12-26-2022 04:04 PM 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^W0(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 RE: [HP41] Lambert function RPN; question - Sylvain Cote - 12-26-2022 04:32 PM (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) RE: [HP41] Lambert function RPN; question - floppy - 12-26-2022 05:46 PM (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 RE: [HP41] Lambert function RPN; question - floppy - 12-26-2022 06:45 PM 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 RE: [HP41] Lambert function RPN; question - Thomas Klemm - 12-27-2022 02:31 AM (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. RE: [HP41] Lambert function RPN; question - Albert Chan - 12-27-2022 04:26 AM (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 RE: [HP41] Lambert function RPN; question - Thomas Klemm - 12-27-2022 04:36 PM 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 RE: [HP41] Lambert function RPN; question - Thomas Klemm - 12-27-2022 04:48 PM (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. RE: [HP41] Lambert function RPN; question - Thomas Klemm - 12-27-2022 05:06 PM (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 RE: [HP41] Lambert function RPN; question - floppy - 12-27-2022 06:42 PM (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 RE: [HP41] Lambert function RPN; question - Pekis - 12-28-2022 11:14 AM 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 RE: [HP41] Lambert function RPN; question - Albert Chan - 12-28-2022 10:26 PM (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^W0 (-0.367879) ≈ 0.368449 321330 e^W-1(-0.367879) ≈ 0.367309 855127 e^W0 (-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 ... RE: [HP41] Lambert function RPN; question - Albert Chan - 12-30-2022 11:27 PM 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^W0 (-0.367879) ≈ 0.368449 321311 e^W-1(-0.367879) ≈ 0.367309 855146 e^W0 (-0.3678) ≈ 0.375551 106237 e^W-1(-0.3678) ≈ 0.360260 737204 RE: [HP41] Lambert function RPN; question - Albert Chan - 12-30-2022 11:48 PM 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