[HP41] Lambert function RPN; question

12262022, 12:13 PM
(This post was last modified: 12262022 06:10 PM by floppy.)
Post: #1




[HP41] Lambert function RPN; question
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/hp41_wlambert/...lambert.41 BUT: What is the line 37 (=RCL M) doing there? I suppose the 2 programs are calculating L(x) in the highbranch 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) HP71B 4TH/ASM/Multimod, HP41CV/X/Y & Nov64d, PILBOX, HPIL 821.62A & 64A & 66A, Deb11 64bPC & PI2 3 4 w/ ILPER, VIDEO80, V41 & EMU71, DM41X 

12262022, 02:07 PM
Post: #2




RE: [HP41] Lambert function RPN; question
(12262022 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. 

12262022, 04:04 PM
Post: #3




RE: [HP41] Lambert function RPN; question
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. (10212020 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) = 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/hp41_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 

12262022, 04:32 PM
(This post was last modified: 12262022 04:32 PM by Sylvain Cote.)
Post: #4




RE: [HP41] Lambert function RPN; question  
12262022, 05:46 PM
(This post was last modified: 12262022 06:11 PM by floppy.)
Post: #5




RE: [HP41] Lambert function RPN; question
(12262022 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 inbetween 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 HP71B 4TH/ASM/Multimod, HP41CV/X/Y & Nov64d, PILBOX, HPIL 821.62A & 64A & 66A, Deb11 64bPC & PI2 3 4 w/ ILPER, VIDEO80, V41 & EMU71, DM41X 

12262022, 06:45 PM
Post: #6




RE: [HP41] Lambert function RPN; question
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 mathprototyping on an HP41 is always fun): Code: 01 LBL "TETR" HP71B 4TH/ASM/Multimod, HP41CV/X/Y & Nov64d, PILBOX, HPIL 821.62A & 64A & 66A, Deb11 64bPC & PI2 3 4 w/ ILPER, VIDEO80, V41 & EMU71, DM41X 

12272022, 02:31 AM
(This post was last modified: 12272022 02:43 AM by Thomas Klemm.)
Post: #7




RE: [HP41] Lambert function RPN; question
(12262022 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 HP42S: Code: 00 { 41Byte Prgm } But it should work for the HP41C 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. 

12272022, 04:26 AM
Post: #8




RE: [HP41] Lambert function RPN; question
(12272022 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=za; 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 

12272022, 04:36 PM
Post: #9




RE: [HP41] Lambert function RPN; question
Nice. This removes another 4 bytes:
Code: 00 { 37Byte Prgm } In addition, it is significantly faster. Example 1E99 XEQ "W" 222.550768956 

12272022, 04:48 PM
Post: #10




RE: [HP41] Lambert function RPN; question
(12262022 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. 

12272022, 05:06 PM
Post: #11




RE: [HP41] Lambert function RPN; question
(12272022 04:26 AM)Albert Chan Wrote: print x Are you seriously still using Python 2? This parrot is dead. Sunsetting Python 2 

12272022, 06:42 PM
Post: #12




RE: [HP41] Lambert function RPN; question
(12272022 04:48 PM)Thomas Klemm Wrote: +1 HP71B 4TH/ASM/Multimod, HP41CV/X/Y & Nov64d, PILBOX, HPIL 821.62A & 64A & 66A, Deb11 64bPC & PI2 3 4 w/ ILPER, VIDEO80, V41 & EMU71, DM41X 

12282022, 11:14 AM
(This post was last modified: 12282022 12:01 PM by Pekis.)
Post: #13




RE: [HP41] Lambert function RPN; question
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 W1 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 

12282022, 10:26 PM
(This post was last modified: 12282022 11:45 PM by Albert Chan.)
Post: #14




RE: [HP41] Lambert function RPN; question
(04102022 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 321330 e^W_{1}(0.367879) ≈ 0.367309 855127 e^W_{0} (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 ... 

12302022, 11:27 PM
Post: #15




RE: [HP41] Lambert function RPN; question
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 321311 e^W_{1}(0.367879) ≈ 0.367309 855146 e^W_{0} (0.3678) ≈ 0.375551 106237 e^W_{1}(0.3678) ≈ 0.360260 737204 

12302022, 11:48 PM
(This post was last modified: 04042023 01:37 PM by Albert Chan.)
Post: #16




RE: [HP41] Lambert function RPN; question
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*(y1/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(1e300, 1)  branch 1 1.4340561272249246e303 697.3227762954601 lua> _ * log(_) 1e300 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 

« Next Oldest  Next Newest »

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