[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: 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: ^{∞}c = W(-ln(c)) / -ln(c) = 1 / e^W(-ln(c))Above quoted example, ^{∞}(√2) = 2We can reuse infinite tetration formula to solve for a^x = b*x Let z = b*x, c = ^{b}√az = a^(z/b) = c^z = c^c^z = ... = ^{∞}chttps://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. 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" 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"` 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 }` 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: 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 }` 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"` 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: +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. 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} (-0.367879) ≈ 0.368449 321330e^W _{-1}(-0.367879) ≈ 0.367309 855127e^W _{0} (-0.3678) ≈ 0.375551 151939 e^W _{-1}(-0.3678) ≈ 0.360260 691185Update: 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^W _{0} (-0.367879) ≈ 0.368449 321311e^W _{-1}(-0.367879) ≈ 0.367309 855146e^W _{0} (-0.3678) ≈ 0.375551 106237e^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)` 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 |