Gamma Function Using Spouge's Method

12182013, 12:55 AM
(This post was last modified: 08052016 06:48 AM by Namir.)
Post: #1




Gamma Function Using Spouge's Method
Gamma Function using Spouge's method.
Memory Map R00 = x and = x1 R01 = a R02 = CHS R03 = Sum R04 = I R05 = sqrt(2*pi) R06 = integer part of I Implementation Code: 1 LBL "GAMMA" 

08072015, 07:23 PM
Post: #2




RE: Gamma Function Using Spouge's Methjod
This is great work, Namir. The Spouge approximation isn't quite as good as the related Lanczos approximation for for the same number of terms in the series, but it is much easier to program, especially on a calculator as the coefficients are simpler.
Les P.S. I have used JMB's GAM+ program in the "old" program listings as my Gamma of choice for years. Will be interesting to see if this performs as well. JeanMichel used the first few terms of standard Stirling series, cleverly factored Hornerstyle to save steps, IIRC. He also stuck to synthetic registers. 

08122015, 04:03 AM
(This post was last modified: 08122015 04:04 AM by lcwright1964.)
Post: #3




RE: Gamma Function Using Spouge's Methjod
Namir,
I looked at the original Spouge paper from 1994 (attached). There is heavy theoretical math in the proofs, and I don't follow it all, but the max relative error is give by a fairly simple formula in Theorem 1.3.1 at the top of page 934. Basically, the larger a is, the smaller the relative error, and the number of terms N = Ceiling(a)  1. You use a =12.5, but a = 13 would do better (in theory) than your choice of 12.5, and you still only do the 12 loops. To compare, Spouge's relative error for a = 12.5 is about 1.2e11, whereas it is about 4.7e12 for a =13. That said it may not make a lot of difference in HP41 or the HP67. I understand that once you start wracking up the terms the promised theoretical precision can get eclipsed by rounding error if one doesn't have a lot of guard digits on hand. I know this is a consideration in arbitrary precision environments (e.g., I read someplace that if you want 45 digits accurate you need to work with 70), but maybe the impact is less with the 13 internal digits of our vintage machines. I don't know if you made your choice of a based on the original paper or another source, but I thought this was interesting. Les 

08122015, 05:39 PM
(This post was last modified: 08122015 06:22 PM by Dieter.)
Post: #4




RE: Gamma Function Using Spouge's Methjod
(08122015 04:03 AM)lcwright1964 Wrote: To compare, Spouge's relative error for a = 12.5 is about 1.2e11, whereas it is about 4.7e12 for a =13. I assume that's the calculated upper bound of the error. But things are different on real calculators. ;) (08122015 04:03 AM)lcwright1964 Wrote: That said it may not make a lot of difference in HP41 or the HP67. I understand that once you start wracking up the terms the promised theoretical precision can get eclipsed by rounding error if one doesn't have a lot of guard digits on hand. That's exactly the crucial point. The theoretical error that would be possible if (!) all calculations were done with infinite (or at least sufficient) precision is much lower than what we can expect from a real calculator with, here, 10 digit working precision. The three additional guard digits you mentioned are only used internally, so they do not matter here. On the '41 or '67, all we got are ten digits. Now let's take a look at a few different values of the variable "a" and see what we get: Mathematically, a value as low as a=9 is sufficient for a 10digit result over the whole (positive) 41/67 working range, i.e. for input between 0 and approx. 71 (overflow limit). Starting at x=0 the error increases with x, reaches a maximum and then decreases again as x gets larger. The largest relative error is 5,33 E–11 near x=32,8. Which is fine for 10 digits (maybe +1 ULP). However, on a real '41 or '67 the last two or even three digits are off. Increasing the value of "a" does not neccessarily help here. For instance, a=13 would be good for a largest relative error of 2,72 E–15 (near x=77,9). But in real life the error gets even larger (!) than for a=9. A value of a=15 theoretically yields an error less than 4 E–17, i.e. full 15–16 digit IEEE double precision. I tried this in Excel VBA as well as on a WP34s in SP mode and got about 10 valid digits near the error maximum. So indeed the working precision has to be greater than the desired accuracy of the result. I assume that 13 digits would suffice for a 10 digit result. But this cannot be accomplished in a simple user program. On the '41, a 13digit MCode implementation would be possible. As long as no additional guard digits are available, a careful choice of "a" is crucial. I would suggest a=9 for the 10digit machines. Larger values do not seem to provide better accuracy. About 8 digits is all we get. Addendum: This also means that a value as high as a=12,5 does not make much sense here. Try Namir's program with a=7 (replace 12.5 with 7 and 1.012 with 1.006). My tests returned results at least as accurate as with a=12,5 (mostly even better) while the program runs about twice as fast. BTW there is an error in the original listing: line 40 should read RCL 06. Dieter 

08122015, 06:03 PM
Post: #5




RE: Gamma Function Using Spouge's Methjod
(08122015 05:39 PM)Dieter Wrote: So indeed the working precision has to be greater than the desired accuracy of the result. I assume that 13 digits would suffice for a 10 digit result. But this cannot be accomplished in a simple user program. On the '41, a 13digit MCode implementation would be possible. Indeed. I believe this is exactly what we have in Angel Martin's SANDMATH. Some years back I was quite engaged with approximating gamma. I learned that Lanczos and Spouge were attractive to some because they were convergent and as such one didn't get into the shift and divide thing for small arguments. Despite this, I think that programmers to this day go with the good oldfashioned Stirling's series, since in practice Spouge and Lanczos don't offer much more in the way of efficiency or accuracy. That said, I really can see that calculator programmers are drawn to the Spouge formula, as the coefficients are computed simply and on the fly. Nonetheless, JM Baillairds excellent GAM+ program in the old HP41 program archive here cleverly rearranges a truncated Stirling approximation in a way that maximizes efficiency and storage of coefficients. Les 

08122015, 09:11 PM
(This post was last modified: 08122015 10:59 PM by Dieter.)
Post: #6




RE: Gamma Function Using Spouge's Methjod
(08122015 05:39 PM)Dieter Wrote: Addendum: This also means that a value as high as a=12,5 does not make much sense here. Try Namir's program with a=7 (replace 12.5 with 7 and 1.012 with 1.006). My tests returned results at least as accurate as with a=12,5 (mostly even better) while the program runs about twice as fast. FTR, here is my implementation for the '41, based on a "hardwired" a=7. Like the factorial function on various other HPs, it actually calculates Γ(x+1). Code: 01 LBL "GAM+1" Since it doesn't seem to get much better than 8 valid digits in the final result, a=7 is sufficient here. Evaluated exactly, the largest relative error is approx. 2,3 E–9. Actually the error may be slightly higher. But it should be less than in the original a=12,5 version. Some examples: Code: 3 XEQ"GAM+1" => 5,999999994 The current version throws an OUT OF RANGE error for x > 54,9. Maybe this can be improved somehow. Edit: the following modification may be a solution. However, the resulting accuracy is slightly less than in the first version. Code: ... This way x=69,9575744 returns 9,999997034 E+99 (exact: ...7560). Dieter 

08192015, 05:32 PM
Post: #7




RE: Gamma Function Using Spouge's Methjod
Thanks, Dieter, this is excellentand very compact too.
Now that I am revisiting the gamma/factorial for calculators issue, it is coming back to me why Spouge and Lanczos, though theoretically attractive because they are convergent and accurate for small real arguments, have seemed less attractive for our calculator work. It is that the just don't do better than the old fashioned asymptotic formulae, even if in theory they should. In Spouge especially, on computes coefficients on the fly and it is no more accurate. I don't think I ever managed to improve on JM Baillard's GAM+ implementation in the old HP41 program archives. His 70step program has some key featuresit uses synthetic registers when storage is needed, does most everything on the stack, cleverly rearranges the Stirling series to simplify the constants embedded in the code, and uses the calculator's FACT function if the argument is an integer. It does the Pochhammer shift for arguments less than 5, which seems to suffice in a 10digit environment to give at least eight good digits and usually nine for most anything. I also rediscovered Peter Luschny's superb page: http://www.luschny.de/math/factorial/app...Cases.html I think I was once drawn to this pseudocode (edited by me to include Luschny's suggested simplification in the formula) which implements the Stieltjes continued fraction: Stieltjes3Factorial(x) y = x+1; p = 1; while y < 7 do p = p*y; y = y+1; od; r = exp(y*(log(y)1)+1/(2*(6*y+1/(5*(y+1/(4*y)))))); if x < 7 then r = r/p fi; return r*sqrt(2*Pi / y) end; However, despite it's compactness and elegance and guarantee of at least 8 digits, I don't think I was ever able to better my favourite Baillard program. Fun times! Les 

08192015, 08:56 PM
(This post was last modified: 08192015 08:57 PM by Ángel Martin.)
Post: #8




RE: Gamma Function Using Spouge's Methjod
(08122015 09:11 PM)Dieter Wrote: This way x=69,9575744 returns 9,999997034 E+99 (exact: ...7560). Yes I used 13digit routines for the MCODE version in the SandMath, which uses the Lanczos formula  and indeed the result for x=70.9575744 is exactly: 9.9999975 99, and using VMANT to see the mantissa: 9.999997560 So it appears to hold up nicely in that region. 

08202015, 11:36 AM
(This post was last modified: 08202015 11:45 AM by Dieter.)
Post: #9




RE: Gamma Function Using Spouge's Methjod
(08192015 08:56 PM)Ángel Martin Wrote: Yes I used 13digit routines for the MCODE version in the SandMath, which uses the Lanczos formula  and indeed the result for x=70.9575744 is exactly: I wouldn't have expected anything less. ;) The accuracy of the Lanczos method, like many others, increases with the argument. The larger x, the more valid digits you get, c.f. Peter Luschny's Gamma/factorial pages. So large values are the easy part. Now the trick is to set a minimum x with sufficient accuracy for which the Lanczos (or another) approximation is applied, and evaluate Gamma for smaller arguments by the usual shift/divide trick. BTW, Sandmath's GAMMA nicely overflows at x=70,95757446 while x=70,95757445 returns 9,999999688 E+99 (exact: ...687). Dieter 

08202015, 12:36 PM
Post: #10




RE: Gamma Function Using Spouge's Methjod
(08192015 05:32 PM)lcwright1964 Wrote: I think I was once drawn to this pseudocode (edited by me to include Luschny's suggested simplification in the formula) which implements the Stieltjes continued fraction: As far as I can see the last x<7 test can be omitted. I tried this in VBA for Excel: Code: Function Gamma(x) The constant in the last line is sqrt(2*pi). It turned out that the results seem to have not just eight but 10 valid digits, at least when evaluated with 15digit precision. Only in some close cases I noticed that the 10th digit may be rounded down instead of up (e.g. for x=1,7 or 7,8). I'll try this one on my '41. ;) Thank you very much for reminding me of Peter Luschny's excellent Gamma/factorial pages. Dieter 

08202015, 07:49 PM
(This post was last modified: 08202015 07:57 PM by lcwright1964.)
Post: #11




RE: Gamma Function Using Spouge's Methjod
(08192015 08:56 PM)Ángel Martin Wrote: Yes I used 13digit routines for the MCODE version in the SandMath, which uses the Lanczos formula May I ask which coefficients? Lanczos offers some in his original paper, and I believe that rendering is quoted in Numerical Recipes verbatim. But there is also the unpublished work of Paul Godfrey, discussed by Viktor Toth on his old http://www.rskey.org website, that reports a matrix algebra method of generating the Lanczos formula coefficients. I think I observed years ago that Godfrey's reported errors in a Matlab implementation were somewhat overly optimistic, as he only seemed to test his formulae with integral and halfintegral arguments. Finally, I think anyone serious about the Lanczos approximation really needs to look at the PhD dissertation of Glen Pugh. Most of the theory and proofs are beyond me, BUT Professor Pugh did some extensive arbitrary precision computer algebra experiments and concluded that the selection of parameters as presented in the original Lanczos paper isn't the best way to go about things, and it is possible to get much less relative error than that originally reported with the same number of terms. For programmers with a mastery of elementary algebra, that part of his work is completely accessible. The Pugh thesis is freely available in a few placessimply Google something like "Pugh Lanczos Dissertation pdf" and you will find it. Dr. Pugh was also receptive to inquirieshe was very friendly with me when I wrote. Hope this is of interest. Les EDIT: Top Google hit for my search: http://web.viu.ca/pughg/phdThesis/phdThesis.pdf Enjoy! 

08202015, 08:22 PM
Post: #12




RE: Gamma Function Using Spouge's Methjod  
08212015, 09:48 AM
Post: #13




RE: Gamma Function Using Spouge's Methjod
(08202015 07:49 PM)lcwright1964 Wrote:(08192015 08:56 PM)Ángel Martin Wrote: Yes I used 13digit routines for the MCODE version in the SandMath, which uses the Lanczos formula Indeed it was the formula with the coefficients published on Viktor's page. You can also see it in the SandMath manual  the most recent version is available at TOS, search for SandMath  and an older version is available here: http://systemyde.com/pdf/SANDMATH44_3x3_Manual.pdf Incidentally I also implemented the same approach on the 41Z module using complex math routines  not 13digit accurate unfortunately but still holding very nicely given the task to perform. It is an ALLMCODE version as well. Thanks for sharing the link to Dr. Pugh's work, very enjoyable. 

08222015, 02:09 AM
Post: #14




RE: Gamma Function Using Spouge's Methjod
(08212015 09:48 AM)Ángel Martin Wrote: Indeed it was the formula with the coefficients published on Viktor's page. Do you mean the seven q's and the formula rearranged for calculators? If that is the case, you might want to know that the selection of the free parameter (Viktor calls it g, but it is a elsewhere) is from the original Lanczos paper and is quite arbitrary. One very practical result of Pugh's work is that there is an optimal value of this free parameter for each number of terms one takes in the series (n or n + 1, depending how you index). Indeed, Pugh published the list as an appendix to his thesis. The optimal value for n = 5 is about 5.58 and is about 6.78 for n = 6. He reports a max relative error of 1.2e10 for n = 5 and 2.7e12 for n = 6, for 9.9 and 11.6 edd (exact decimal digits) respectively. Keep in mind that these errors cover complex arguments as well, and the performance for eligible real values tends to be somewhat better. I can swear that I did up a Maple worksheet of this ages ago, that accepted these optimal parameters as input, solved Godfrey's matrix equations for the associated coefficients, doing it all to very generous arbitrary precision, and transformed them into the form recommended by Victor for calculator programming. But can I find the darn thing? Alas, nothat's a couple of hard drives ago. If I do redo this I can compute and send you "improved" coefficients, in the event you ever tweak the code. It might not make much difference unless you can shorten your series by a term, but I sure would like to ponder the question. I want to continue this discussion elsewhere in the forum as I have some broader newbie questions about the newest SandMath, Lib#4, and programming Nov64 with this stuff (if possible). Les 

08222015, 08:35 AM
Post: #15




RE: Gamma Function Using Spouge's Methjod
(08222015 02:09 AM)lcwright1964 Wrote:(08212015 09:48 AM)Ángel Martin Wrote: Indeed it was the formula with the coefficients published on Viktor's page. Well, yes I tweak the code all the time (in fact just recently since Dieter found a bug in FLOOR)  but first of all, the reference I use doesn't have any free parameter that I can see. This is the link, and the coefficients are already shown in the more efficient way, I think: http://www.rskey.org/CMS/index.php/thelibrary/11 I'm sure I'm missing something so will appreciate the enlightenment. SandMath questions are always welcome, as a matter of fact I'm goinf to burn my own NoVRAM today with the updated version. It's been a while I did this, I mostly work on my CL units of course  much more comfortable and super fast. \AM 

08222015, 08:58 AM
Post: #16




RE: Gamma Function Using Spouge's Methjod
(08222015 08:35 AM)Ángel Martin Wrote: Well, yes I tweak the code all the time (in fact just recently since Dieter found a bug in FLOOR)  but first of all, the reference I use doesn't have any free parameter that I can see. This is the link, and the coefficients are already shown in the more efficient way, I think: The free, or arbitrary, parameter used in the Numerical Recipes Lanczos formula and rearranged by Viktor is embedded in the 5.5. It is little gamma in the original paper, g on Viktor's page, and r in the Pugh thesis. In the formula you are using it happens to be 5. Read further down on that very page, to the discuss of the Godfrey matrix decomposition, where Viktor identifies this g as "an arbitrary real parameter." The key to minimizing relative error to choose this parameterand yes it is chosen, and the coefficients then computed in the complex way according to it and the number of desired termsjust so. Both Lanczos and Godfrey, who inspired Viktor's calculations, were more or less arbitrary in the selection. Pugh did empirical tests and theoretical analysis to find the best value for a desired number of terms. Indeed, in subsequent years he extended this work up to n = 100 (the thesis only goes to n = 60) and I just found that page that he sent me in 2006. In it he gives the parameters to more digits and has improved the error estimates by using a more rigourous calculation. And indeed, it seems that in the very situation you use (n = 6, or 7 coefficients with the index starting at 0), the optimal value of g is around 6.78, not 5. That's what I mean by "free parameter." Sorry for the confusion. As for SandMath use, I am curious about two things: 1. What is Library 4 and why do I need it? 2. Given the switching going on in SandMath (as in the Advantage Pac), how do I correctly configure ClonixConfig to burn it into my Nov64? Thanks again, Les 

08222015, 09:57 AM
(This post was last modified: 08242015 08:42 AM by lcwright1964.)
Post: #17




RE: Gamma Function Using Spouge's Methjod
Well I did find that Maple worksheet. Looks like I basically ported Viktor's C++ code for Godfrey's method and took advantage of Maple's easy to use linear algebra capabilities.
Unfortunately I did NOT expand things out and collect terms as Viktor did (and you used) so my result looks a bit different. My preliminary conclusion is that going with Pugh's recommendation for the arbitrary parameter one can keep the relative error the same (at least 10 digits good for positive real arguments as long as one computes in 13digit precision), with one fewer terms in the series. Here is my result with the coefficients presented to 13 significant digits: (0.9446965718570e2 + 1.367185253712/z  1.839535625335/(z+1) + .6838535775175/(z+2)  0.6542905498252e1/(z+3) + 0.6572071053721e3/(z+4)) * ((z+5.081000210223)/exp(1))^(z1/2) Input has been adjust to compute Gamma(z), not Gamma(z+1) or z!. The scaling by e in the second factor before exponentiation should prevent overflow when other formulas tend to for larger arguments, but I have not rigourously examined that. In this case the arbitrary parameter recommended by Pugh is 5.58100021022. In the original paper Lanczos does not present an example with n = 5 (so 6 coefficients) for comparison. In the version you use now you embed 7 12digit coefficients and one 2digit parameter5.5 (from 5 + .5). Here you embed 6 13digit coefficients and one 13digit parameter. And there is one less series term. So this should be faster and smaller, at least a little. I should work on doing what Viktor does so as to produce a polynomial divided by a product, etc. In the meantime, this is what I have. I work in Mathematica more these days so I should see what happens there. Les 

08222015, 10:45 AM
(This post was last modified: 08222015 10:49 AM by Ángel Martin.)
Post: #18




RE: Gamma Function Using Spouge's Methjod
(08222015 08:58 AM)lcwright1964 Wrote: As for SandMath use, I am curious about two things: Thanks for the additional explanation on the parameter's role in the formulas, I think I can follow now. The Library#4 is a ROM that contains a library of subroutines used by several other modules of mine. This saved lots of space and made the code/operation much more consistent across all modules. It sits in page#4, where typically there's nothing plugged  exceptions being the diagnostics ROM and a disabled printer on the HPIL module. The modules that use it are as follows: (I think I'm not missing anything but it's from memory)  OSX/3 and OSX/4  one page w/ 4 banks in it = 4 ROMs CLUTILS and POWERCL  one page w/ 4 banks in it = 4 ROMs ALPHA44 RAMPAGE 4L TOOLBOX 4 41Z Complex Numbers  2 pages = 2 ROMS SANDMATH  2 pages of 4 banks each = 8 ROMS SANDMatrix  2 pages of 2 banks each = 4 ROMs TotalREKALL HP16C Emulator (about to be released)  one page w/ 4 banks = 4 ROMs ClonixConfig has a very flexible U/I that allows you to burn bankswitched modules like the OS/X, SandMath & SandMatrix. You simply need to select the right combination of banks/pages more or less in contiguous sequence. The attached pic shows an example of the SandMath + SandMatrix burned on a ClonixD with JOINED blocks. A Novo64 would use the same setup, once you disable the HEPAX RAM checkbox. This doesn't get you there yet  as it's still missing the Library#4. So you'd probably need another Clonix41 with the Lib#4, and (highly recommended) the OS/X4. You can throw in there the TotalRekall to fill it up to 6 ROMS. Hope this clarifies, let me know if you have more questions. 

08222015, 01:17 PM
(This post was last modified: 08222015 02:31 PM by Ángel Martin.)
Post: #19




RE: Gamma Function Using Spouge's Methjod
forgot to mention that all those modules (Lib#4 included) are available at TOS so you can do your dry runs on V41 to get a pretty good idea of how things work and what's in there.
ps. see two suggestions for Clonix41 configurations as per previous post 

08222015, 11:51 PM
Post: #20




RE: Gamma Function Using Spouge's Methjod
Thanks for all this. I have to admit I have obviously been out of touch. The last time I paid attention to SandMath I think it consisted of just lower and upper halves, and that's about itbasically, what is still included with Diego's most recent ClonixConfig software. I took a peak at the most recent offerings on TOS, including the simpler SandMathIV module that doesn't use Library #4, and modfile.exe tells me that there is a lot more going on with these things in the three or four years since.
And thanks for the configuration suggestion. I would probably just burn the SandMath 4x4 constituents with Library #4, and leave SandMatrix and the others be for now. Of course, you have tempted me to order another of the Clonix family from Diego Getting back to the original topic, I have used that very same Maple worksheet and Pugh's recommendations to give this n = 4 (i.e. 5 coefficient) version for a 10digit FOCAL environment: (0.3264892413e1 + 1.188688981/z  1.068410309/(z+1) + .1887119902/(z+2)  0.2745021709e2/(z+3)) * ((z+ 3.840881909)/exp(1))^(z1/2) In an ideal situation (i.e., a couple of guard digits kicking around) this is supposed to give between 9.5 and 11 edd in the z = 0 to 70 range, but of course we don't have guard digits so 8 to 10 is about the best I think it would do. Given this I think JMB's implementation of the Stirling series in GAM+, with its simpler coefficients, will still do the same or better. That said, I would like to see how this fares in FOCAL code for the HP41. It would even work for the 67/97 where one could store the program on one side of the card (it would be very short) and the six 10digit constants as register data on the other side . Without Pochhammer shifting for small arguments this should be faster than a Stirling approach, and of course the Spouge approach uses more terms and calculates more. If I actually turn this into program steps for the 41 or the 67/97 I will post. Les P.S. I notice that the coefficients (q0 to q6) from Viktor's page are only presented to 12 digits, not 13. Is that what you are using, or did he or you recompute 13digit versions for use in your MCODE? That does seem possible given that the original Lanczos coefficients, p0 to p6, are given just above them to up to 16 digits. 

« Next Oldest  Next Newest »

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