From comp.sys.handhelds Tue Feb 5 09:10:38 1991 Path: mentor.cc.purdue.edu!purdue!bu.edu!rpi!zaphod.mps.ohio-state.edu!swrinde!cs.utexas.edu!uunet!shelby!agate!dog.ee.lbl.gov!usenet From: austern@ux5.lbl.gov (Matt Austern) Newsgroups: comp.sys.handhelds Subject: Gamma function Summary: Here's a program to compute gamma(z) for arbitrary complex z. Message-ID: <9483@dog.ee.lbl.gov> Date: 4 Feb 91 06:17:59 GMT Reply-To: austern@ux5.lbl.gov (Matt Austern) Organization: Lawrence Berkeley Laboratory (theoretical physics group) Lines: 51 The gamma function is one of the standard "special functions" of complex analysis. It is therefore a little bit frustrating to see that the gamma function on the HP 48 only accepts real arguments. (Not that I'm complaining, of course: HP is the only company I know of that puts the gamma function on its calculators at all.) The following program corrects this. The algorithm isn't as accurate as that on the HP 48, but it still is quite good: better than 1 part in 10^9. This program uses an algorithm from Numerical Recipes. (Press, et. al.) They give no theory behind the algorithm, besides the obvious observation that it asymptotically approaches Stirling's approximation, but they do give a reference. It is only valid for Re (z) > 1, so I use the "reflection formula" if given an argument with Re(z) < 1. You'll note that this program calls itself (in the reflection formula), so you should either store it under the name GAMMA or change that call. %%HP: T(3)A(R)F(.); \<< \-> z1 \<< z1 1 - \-> z \<< IF z RE 0 < THEN \pi \->NUM z * DUP SIN / 1 z - GAMMA / ELSE z 5.5 + z .5 + ^ z 5.5 + NEG EXP * 2 \pi \->NUM * \v/ * 1 76.18009173 z 1 + / + -86.50532033 z 2 + / + 24.01409822 z 3 + / + -1.231739516 z 4 + / + .00120858003 z 5 + / + -.00000536382 z 6 + / + * END \>> \>> \>> -- Matthew Austern austern@lbl.bitnet Proverbs for paranoids, 3: If (415) 644-2618 austern@lbl.gov they can get you asking the wrong questions, they don't have to worry about answers.