SandMath routine of the week: Inverse Gamma Function  Printable Version + HP Forums (https://archived.hpcalc.org/museumforum) + Forum: HP Museum Forums (https://archived.hpcalc.org/museumforum/forum1.html) + Forum: Old HP Forum Archives (https://archived.hpcalc.org/museumforum/forum2.html) + Thread: SandMath routine of the week: Inverse Gamma Function (/thread240692.html) 
SandMath routine of the week: Inverse Gamma Function  Ángel Martin  03152013 As a belated celebration of piday, here's a short routine to calculate the Inverse Gamma function (positive arguments only)  Simply using Newton's method and the SandMath, it doesn't get any shorter  and on the CL it quickly converges after a few iterations... worth seeing it in action!
01 LBL "IGMMA" Example: in addition to x=3, which other argument also returns Gamma(x) = 2 ?
Answer: 2, XEQ "IGMMA" > 0,4428773962
Cheers, Note: it's a rather crude firstpass, go ahead and improve on it if you're interested...
Edited: 15 Mar 2013, 10:48 a.m.
Re: SandMath routine of the week: Inverse Gamma Function  Walter B  03152013 Buenas tardes, Ángel, what do you use PSI for in your function set? TIA for enlightenment.
d:?
Re: SandMath routine of the week: Inverse Gamma Function  Ángel Martin  03152013 The above example shows one usage. PSI is a handy trick to obtain the derivative of Gamma, as well as being very useful in probability problems.
Re: SandMath routine of the week: Inverse Gamma Function  Gerson W. Barbosa  03152013 Very nice short routine! %%HP: T(3)A(D)F(.); Cheers, Gerson.
Re: SandMath routine of the week: Inverse Gamma Function  Paul Dale  03152013 We should have included PSI in the 34S too...
Re: SandMath routine of the week: Inverse Gamma Function  Ángel Martin  03162013 Glad it piqued your curiosity, Gerson. Interestingly enough the SandMath implementation does converge for x=3 and x=4, albeit it takes a long time. It's somewhat of a miracle because the estimations become negative for a while and we all know what happens to Gamma in the negative plane. I guess there must be a clever way to tweak the calculation to reduce the number of iterations, perhaps keeping the action in the positive plane. In a google search I found almost no references to the InverseGamma function, so I guess it's totally uninteresting (which makes it attractive to my eyes :) Best,'AM
Edited to include the links:
Edited: 16 Mar 2013, 6:21 a.m.
Re: SandMath routine of the week: Inverse Gamma Function  Walter B  03162013 Something left for the 43S ;) Anyway it's included as a library function AFAICS.
d:)
Re: SandMath routine of the week: Inverse Gamma Function  Paul Dale  03162013 Yes, it is in the standard library. It was also coded as an internal function but it got left out :( The 43S will definitely have this one. A native inverse gamma is tempting too :)
 Pauli
Re: SandMath routine of the week: Inverse Gamma Function  Gerson W. Barbosa  03162013 ¡Hola, Ángel!
Quote: References to InverseGamma are indeed hard to find. I think I had already stumbled onto your first link last year when I was trying to implement this function. The best I was able to to gave only three digits. The method you have presented, using Psi, is the best I've ever seen. Thanks again for posting! Here is a shorter and faster RPL version. The number of iterations is 6 or 7 for most arguments. It will work in the range 1 < x <= 1.34830732701E488. Best regards,
Gerson. InvGamma: Re: SandMath routine of the week: Inverse Gamma Function  Ángel Martin  03162013 Hi Gerson, I really have trouble to read RPL, can you provide the RPN equivalent for the calculation of the initial guess? Looks like a substantial improvement to the initial idea! This is what I figured out:
01 LBL "IGMMA" Best, 'AM
Edited: 16 Mar 2013, 6:39 p.m.
Re: SandMath routine of the week: Inverse Gamma Function  Gerson W. Barbosa  03162013 Hi Ángel, OVER is equivalent to RCL Y, but it's easier if you know where these constants and expression come from (see below *).
This can be a possible RPN equivalent for the first part of the program: 01 LBL "IGMMA" Cheers, Gerson. 
(*) The following is part of a curve fitting made with help of an old version of
_{ Model Definition: Y = x/(a+b*x+c*sqr(x)) Number of observations = 172 Number of missing observations = 0 Solver type: Nonlinear Nonlinear iteration limit = 250 Diverging nonlinear iteration limit =10 Number of nonlinear iterations performed = 23 Residual tolerance = 0,0000000001 Sum of Residuals = 1,97982012453611E02 Average Residual = 1,1510582119396E04 Residual Sum of Squares (Absolute) = 0,071461272245453 Residual Sum of Squares (Relative) = 0,071461272245453 Standard Error of the Estimate = 2,05632625029686E02 Coefficient of Multiple Determination (R^2) = 0,9994754674 Proportion of Variance Explained = 99,94754674% Adjusted coefficient of multiple determination (Ra^2) = 0,99946926 DurbinWatson statistic = 0,181006272327966 Regression Variable Results Variable Value Standard Error tratio Prob(t) a 0,194180061078165 6,75779102290044E02 2,873425065 0,00458 b 0,160140574611016 5,79167114529964E04 276,5014977 0,0 c 2,20860298743044 1,42205819324635E02 155,3103099 0,0  X Value Y Value Calc Y Residual % Error Abs Residual Min Residual Max Residual 1 0 0 0 0 0 0 0,09225268066 0,09044327801 2 0,6931471804 0,2310490602 0,3233017409 0,09225268066 39,92774547 0,09225268066 3 1,791759469 0,4479398673 0,5212429458 0,07330307854 16,36449084 0,07330307854 4 3,17805383 0,635610766 0,6848643811 0,0492536151 7,749021529 0,0492536151 5 4,787491743 0,7979152904 0,8263771459 0,02846185547 3,567027204 0,02846185547 6 6,579251212 0,9398930303 0,9517404203 0,01184739003 1,260504084 0,01184739003 7 8,525161361 1,06564517 1,064573069 0,001072101061 0,10060582 0,001072101061 ... 167 691,1834011 4,114186911 4,091157824 0,02302908714 0,5597481991 0,02302908714 168 696,3073651 4,120161924 4,096382192 0,02377973199 0,5771552776 0,02377973199 169 701,4372638 4,126101552 4,101568196 0,02453335547 0,5945892306 0,02453335547 170 706,5730622 4,132006212 4,106716318 0,02528989421 0,6120487945 0,02528989421 171 857,9336698 4,289668349 4,241416568 0,04825178142 1,124837109 0,04825178142 172 1128,523771 4,514095084 4,423651805 0,09044327801 2,003574943 0,09044327801  } In the X column, ln(InvGamma(x)); in the Y column, ln(InvGamma(x))/x. This was the best fit out of 57 models.
Edited: 16 Mar 2013, 9:52 p.m.
Re: SandMath routine of the week: Inverse Gamma Function  Ángel Martin  03172013 Hi Gerson, many thanks for the RPL help and the additional details on the data fit  yes, that's a good idea that has prompted me to look for other methods to get the initial estimation, with the goal to reduce the number of iterations even further (albeit your method is already very good, perfect for fast machines like the 50G and the CL :) I thought about some value based on the Stirling approximation, and that took me straight to the Lambert function. Now I can see how things fit in the first URL, so I'm using David Cantrell's formula as initial guess. This works for x greater than x0=1.461632, the zero of PSI (aka the local minimum of Gamma for x>0). I've added a branch to deal with x<x0, also including arguments up to x=2  somehow I like the "alternative" result more than the standard x=3. This branch goes to the negative values for the solution and takes a lot of iterations, but still it's better that nothing. Below I've listed both versions, the one on the right side is using your approach. The other is two lines longer but uses the same room in memory. Indeed the number of iterations is reduced, saving precious time in slow machines  despite the initial frontloaded calculation of WL0. Somehow is nice to have these three functions (Lambert, PSI and Gamma) in the same code, a nice trinity working together. Needless to say this will make it in the next revision of the SandMath, which BTW will have two bankswitched banks for each page  four in total.
01 LBL "IGMMA" LBL "IGMMA"
Cheers,
Edited: 17 Mar 2013, 6:16 a.m.
Re: SandMath routine of the week: Inverse Gamma Function  Gerson W. Barbosa  03172013 Hi Ángel,
Quote: Wouldn't it be better to write s specific program suited for that branch? Perhaps only the initial estimates should be different. It is interesting to notice that PSI can be replaced with LN for the positive branch (arguments > 1, returning values > 2), since Psi(x) get closer to ln(x) as x grows. Thus, an HP42S version is possible. There will be occasional differences of one unit it the last significant digits for arguments < 30 or so. There is a couple of unnecessary operations in the RPL program above (also in the RPN version). Fix below. Cheers, Gerson. 
InvGamma:
01 LBL "IGMMA" P.S.: This HP42S version uses only the stack. There are no noticeable differences in the running times however. The conversion to the wp34s should be straightforward.
00 {74Byte Prgm } 19 Rv 
Edited to remove my comment about your use of Lambert W. Somehow I imagined it was inside of loop. I think I hadn't woken up yet :)
Edited: 18 Mar 2013, 12:52 a.m.
Re: SandMath routine of the week: Inverse Gamma Function  Ángel Martin  03182013 Very nice Gerson, I think we've nailed it  save the "left" branch of course.To be continued... :) I went ahead and timed the execution times for both versions, see the table below. Using the more complex initial estimation yields consistently shorter times as you can see. Times in seconds.Never mind the actual values (on V41, standard speed) but the deltas tell the story:
x D. Contrell DataFit
Cheers, Re: SandMath routine of the week: Inverse Gamma Function (WP 34S version)  Gerson W. Barbosa  03182013 Hi Ángel,
Significantly faster! Better use the curve fit equation only in programs for calculators that lack the Lambert W then. Cheers, Gerson.
WP 34S version: Re: SandMath routine of the week: Inverse Gamma Function (WP 34S version)  Paul Dale  03182013 Quote:
XEQ'[PSI]' This can be entered via the keyboard and the gshift Greek letters or via CAT and navigating there and pressing XEQ.  Pauli
Edited: 19 Mar 2013, 2:27 a.m.
Re: SandMath routine of the week: Inverse Gamma Function (WP 34S version)  Paul Dale  03182013 Nice use of some of the 34S's new features by the way. Can I add this routine to the library?
Re: SandMath routine of the week: Inverse Gamma Function (WP 34S version)  Gerson W. Barbosa  03182013 Of course! Again, thanks Ángel for presenting us the method :)
Gerson.
Re: SandMath routine of the week: Inverse Gamma Function (WP 34S version)  Ángel Martin  03192013 Gerson you're a gentleman  my pleasure of course, and thanks to you for the improvements.
Re: SandMath routine of the week: Inverse Gamma Function  peacecalc  03212013 Hello all, nevertheless it is a interesting algorithm, but the hp50g can solve the problem with it's own built in function "ROOT", a very nice and short RPL proggie:
Interesting feature: If you use a 1 as a initial guess, you get the negative solutions.
Greetings
Re: SandMath routine of the week: Inverse Gamma Function (WP 34S version)  Gerson W. Barbosa  03222013 Compatible with SSIZE8 and one step shorter:
+ LBL'[GAMMA][^1]' I guess [<>] XZZY takes longer than y[<>] T, but I haven't noticed any difference when the argument is 0.89 (21 second with a chronometer, about the same time I obtained with TICK in the previous version). DBLON would require the CVNG?02, but I don't see a way to do this without extra instructions inside to loop. Simply changing the parameter to 02 would cause endless loop for certain arguments.
Gerson.
Re: SandMath routine of the week: Inverse Gamma Function  Gerson W. Barbosa  03222013 Hi,
Interesting technique for using the numeric solver in a program. Cheers,
Gerson.
Re: SandMath routine of the week: Inverse Gamma Function (WP 34S version)  Paul Dale  03222013 Gerson, thanks for this update. I suspect there isn't much time difference between the shuffle command and the swap. The other overheads of interpreting instructions will likely far exceed the time to do either. What about using CNVG? 03. This picks CNVG? 00 in single precision mode and CNVG? 02 in double.
Re: SandMath routine of the week: Inverse Gamma Function (WP 34S version)  Gerson W. Barbosa  03232013 Thanks! I'd found this information in page 91 of the Blue Book:
"3 = Choose the best for the mode set, resulting in taking 0 for SP and 2 for DP (see Appendix H)." This is plain clear, but the reference to Appendix H somehow diverted me. I was expecting to find more information about CNVG?, but Walter was meaning more information about SP and DP, I realize now. So, let it be CNVG? 03 then. Gerson.
Re: SandMath routine of the week: Inverse Gamma Function  peacecalc  03242013 Hello Gershon,
Original from Gerson W. Barbosa: ... a specific solver will usually be faster Of course that is right, but it is easy to use (without a great bunch of knowledge in numeric analysis). And I've stated an interesting behavior of my "brutal force" method: With huge arguments my 50g needs let's say some seconds for the solution, and with small arguments in the scale of 10, it is very fast compared to the values of post #15 that posted Gershon. Certainley the duration of computation depends on the displaying mode because I use the built in function "ROOT".
Greetings Re: SandMath routine of the week: Inverse Gamma Function (WP 34S version)  Gerson W. Barbosa  03242013 Thanks! As I suspected, the use of LN instead of Psi is responsible for the overall slowness, especially for arguments in the beginning of the range. Calling the library digamma program would take up only one step, but it doesn't preserve the stack. This would have to be taken care of, which would require extra steps. Since Psi need not be exact in this application, I decided to use an approximation:
Psi(x) ~ ln(x)  1/(2*x)
Edited: 25 Mar 2013, 12:20 a.m.
Re: SandMath routine of the week: Inverse Gamma Function (WP 34S version)  Paul Dale  03252013 Quote: The diacritic is gone now.
Quote: In a library routine, I think so. If anyone objects, they can always key in the earlier version or delete the routine completely.
Re: SandMath routine of the week: Inverse Gamma Function  Gerson W. Barbosa  03252013 Here is another HP42S version. Psi(x) is approximated as ln(x)  1/(2*x) instead of simply ln(x) to speed execution up, especially for small arguments.
00 {86Byte Prgm } 23 STO ST Z Re: SandMath routine of the week: Inverse Gamma Function (WP 34S version)  Gerson W. Barbosa  03282013 Saving two instructions now, so five extra steps actually:
LBL'[GAMMA][^1]' Gerson.
Edited: 28 Mar 2013, 7:19 p.m.
Re: SandMath routine of the week: Inverse Gamma Function (WP 34S version)  Gene Wright  03282013 Fabulous work, Gerson!
Re: SandMath routine of the week: Inverse Gamma Function (WP 34S version)  Gerson W. Barbosa  03292013 Not quite! I completely forgot to test it in SSIZE8 mode. When I finally did it this afternoon, I noticed it was taking much longer than when in SSIZE4 mode. Fortunately the fix is very easy: just replace the only RCL T instruction with RCL Z:
... Occasional stacksizewise incompatibilities are a side effect of having two stack sizes. I only change to SSIZE8 when doing many calculations involving complex number (which lately I have no need to). I would prefer a parallel fourlevel stack for doing this, though, as in the HP15C. But there's a reason this couldn't be implemented on the wp34s, so the 8level stack for such situations is a blessing (but only in this case, IMHO). I wonder what would be the consequences of usersettable stack size, as some have suggested... A really fabulous work was yours, adapting matrices routines to work on the wp34s! Regards,
Gerson.
Re: SandMath routine of the week: Inverse Gamma Function (WP 34S version)  Walter B  03292013 Obrigado, Gerson, por a sua programa! Quote:All a matter of habitude: I run on SSIZE8 all the time. I love RPN calculating free of stackoverflow worries.
d:)
Re: SandMath routine of the week: Inverse Gamma Function (WP 34S version)  Gerson W. Barbosa  03302013 Quote:
De nada!
Quote: Yes, I know there is a number of wp34s user who prefer the 8level stack. That's why I realized all library programs should ideally run under both stack sizes, without changing the user setting. It would be easier just to save the user settings in the beginning of the program, set the proper stack size and restore the user settings when the program returns, but I don't know whether this is possible. At least I haven't found a way to do it (but I haven't read all the book yet). Regards,
Gerson.
Re: SandMath routine of the week: Inverse Gamma Function (WP 34S version)  Paul Dale  03302013 Making a program behave like a built in command is fairly easy. It is much easier in XROM but that isn't available to users.
LBL'ABC' And despite Walter's comment an eight level stack can still run out, especially when programming :)  Pauli
Edited: 30 Mar 2013, 2:35 a.m.
Re: SandMath routine of the week: Inverse Gamma Function (WP 34S version)  Walter B  03302013 Pauli, Quote:If you want to sabotage a system you can always, of course ;) As stated above, I was talking about calculating not programming. Some minimum intelligence given, you cannot overload an eightlevel stack in real calculations IMHO.
d:)
Re: SandMath routine of the week: Inverse Gamma Function (WP 34S version)  Gerson W. Barbosa  03302013 Thanks! I would use this technique in complex and long programs. Preserving the stack only is a nice bonus. However, in a relatively short program like Gamma^{1} seven extra steps might be too many.
Quote: I don't remember ever getting out of stack when doing calculations, but then again mine usually are not complicated. The 8level stack should be hand for calculating this one, though:)
Gerson.
Re: SandMath routine of the week: Inverse Gamma Function (WP 34S version)  Walter B  03302013 Quote:Hmmmh  let me see: input x y z tUnless I'm mistaken, three levels are sufficient here.
d:?
Re: SandMath routine of the week: Inverse Gamma Function (WP 34S version)  Gerson W. Barbosa  03302013 Quote:
Indeed. Your result is correct: 0.835724531752514I remember there is an expression somewhere in the archives which is not solvable within four stack registers, unless intermediate results are stored in a numbered register, but I cannot find it right now. This example is from one of many discussions about RPN versus AOS:
Mach number  RPN vs. Algebraic Edited: 30 Mar 2013, 6:16 p.m.
Re: SandMath routine of the week: Inverse Gamma Function (WP 34S version)  Thomas Klemm  03302013 Quote: Maybe this? Limitations of RPN (was: HP32s Engineering Applications Manual) But there was an earlier thread: It Can Be Done!
Kind regards Re: SandMath routine of the week: Inverse Gamma Function (WP 34S version)  Gerson W. Barbosa  03312013 Thanks, Thomas! Yes, that was the expression you presented in message #14 in the first link. Of course no problem with it on the WP34S with SSIZE8 mode on. Your comment is very impressive: "Whichever solution you might choose you will always run into a stackoverflow and what's worse you probably won't even notice." Walter does have a point :) Cheers, Gerson.
