|Enhanced version optimized for 12cp|
Message #5 Posted by Les Wright on 31 Jan 2009, 4:25 a.m.,
in response to message #1 by Les Wright
Here is my enhanced routine, where the initial 12 steps carry out the the Pochhammer shift-and-divide correction. Basically, the argument is incremented until it exceeds the integer threshold stored in register 2 (I have been experimenting with 7 thru 10, and seem to like 9 so far), computes the gamma of the incremented argument, and divides out the accumulated Pochhammer product in register 3 at the end. If anyone is familar with Jean-Marc Baillard's routine GAM+ in the HP41 library, I have simply adapted his approach to the Pochhammer shift to the limitations of the 12CP, but, unlike him, I have stuck with the Stieltjes continued fraction for the main computation (whereas he uses a traditional Stirling series):
The usage is similar to that reported in my initial post. I store the 12-digit version of sqrt(2Pi) in register 1 and the desired threshold in register 2 (in my case, this means 9 STO 2). And then I compute gamma for positive values exactly as described earlier. Some examples, once I have stored sqrt(2Pi) in R1, 9 in R2, and reset the program pointer with f CLEAR PRGM (in run mode):
.25 R/S returns 3.625609908 (all 10 digits correct)
5 R/S returns 24.00000000 (all correct)
3.222 R/S returns 2.478039436 (actual is 2.478039437)
5.5 R/S returns 52.34277778 (all correct)
.5 R/S returns 1.772453851 (this is sqrt(Pi), all digits correct)
.0001 R/S returns 9999.422884 (actual is 9999.422883)
12.771 R/S returns 269174672.9 (actual is 269174673.0)
Speed is really good too, though the smaller arguments take a little longer since the Pochhammer shift loop takes a few more trips in its augmentation of the argument. Of course, larger arguments that don't need shift lead to quicker results.
The routine WILL handle negative arguments since the Pochhammer shift steps will move the argument into positive territory, but for negative arguments with larger absolute value this can be rather inefficient. For example, -65.126 R/S will return the correct result of 5.831779435E-91, but it is much slower than for positive values since the Pochhammer shift loop has to go around seventyish some odd times to get the argument above 9. (That said, we are talking only a dozen or so seconds as opposed to a couple.) Programming the sine function to use the reflection formula, as Egan did, is much more efficient.
Also, keep in mind that the impressive accuracy is an advantage of the 12CP over the 12C--stack and register computations use 12-digit numbers all the way through, even though only 10 are displayed. A 12C version of this would NOT give as impressive results since rounding errors affect visible digits, not "hidden" ones.
Folks who insist that the routine return the factorial just like the f x! key sequence on the 15C need only increase the argument by 1, either manually or right at the beginning of the program, since of course Gamma[x+1] = x!
I really love programming the oft maligned 12cp, since I find the 12-digit accuracy and speed of the calculator are strong advantages (though the stubborn keyboard is not). I hope folks enjoy this and will see it as a supplement or alternative to Egan's excellent contribution.
Edited: 31 Jan 2009, 4:44 a.m.