Simon Jeffery's Software Store

TAP: mathematical functions.

Quick Reference

FUNCTION TAP_GAMMA ( X, STATUS )	Evaluate the gamma function Gamma (x) for a real argument.
FUNCTION TAP_GAMMLN ( X, STATUS )	Evaluate the gamma function ln Gamma (x) for a positive real argument.
FUNCTION TAP_FACTRL ( N , STATUS )	Return the value of factorial N, i.e. N! as a floating point number
FUNCTION TAP1_POWERG ( Z, STATUS )	Evaluate power series expansion for ln Gamma (z)

---

Long Reference

---

DOUBLE PRECISION FUNCTION TAP_GAMMA ( X, STATUS )

Purpose:
   Evaluate the gamma function Gamma (x) for a real argument.

Invocation:
   DOUBLE PRECISION TAP_GAMMA
   RESULT = TAP_GAMMA( X, STATUS )

Description:
   The gamma function (Gamma (x)) is evaluated for a single real
   argument.

   The gamma function for positive x using the TAP function
   TAP_GAMMLN, which uses a power series expansion. If x < 0
   the reflection formula

       Gamma(x) * Gamma(1-x) = pi.cosec(pi.x)

   is applied. If x = -n, n = 0,1,2,..., Gamma(x) has a pole.
   An error status is returned, and TAP_GAMMA (x) is not defined.

Arguments:
   X = DOUBLE PRECISION (Given)
      The real number X for which GAMMA(X) is required.
   STATUS = INTEGER (Given and Returned)
      The global status.

Returned Value:
   TAP_GAMMA = DOUBLE PRECISION
      GAMMA(X) as a real number

Subprograms used:
   DOUBLE PRECISION FUNCTION TAP_GAMMLN ( Z, STATUS )
   DOUBLE PRECISION FUNCTION TAP1_POWERG ( Z, STATUS )

---

DOUBLE PRECISION FUNCTION TAP_GAMMLN ( X, STATUS )

Purpose:
   Evaluate the gamma function ln Gamma (x) for a positive real argument.

Invocation:
   DOUBLE PRECISION TAP_GAMMLN
   RESULT = TAP_GAMMLN( X, STATUS )

Description:
   The gamma function (Gamma (x)) is evaluated for a single real
   positive argument, the logarithm (ln (Gamma(x))) is returned.

   The gamma function is evaluated using a power series expansion
   valid over the interval 0 < x < 2. However, the expnasion is only
   used in the interval 1/2 < x < 3/2, where it is most rapidly
   convergent. For other values of the argument x, a combination of 
   the reflection formula

       Gamma(x) * Gamma(1-x) = pi.cosec(pi.x)

   and the recurrence relation

       Gamma(x+1) = x.Gamma(x)

   are applied. The algorithm is not as fast as obtained with
   a series approximation, (eg Stirling's formula or the Lanczos
   approximation - see Press et al.), but is adopted here
   in order to achieve maximum arithmetic precision.

References.
   Abramowitz M., &  Stegun I.A., 1974, Handbook of Mathematical Functions,
      Dover.
   Press et al. 1989. Numerical Recipes  

Arguments:
   X = DOUBLE PRECISION (Given)
      The real number X for which GAMMA(X) is required.
   STATUS = INTEGER (Given and Returned)
      The global status.

Returned Value:
   TAP_GAMMLN = DOUBLE PRECISION
      LOG (GAMMA(X)) as a real number

Subprograms used:
   DOUBLE PRECISION FUNCTION TAP1_POWERG ( Z, STATUS )

---

DOUBLE PRECISION FUNCTION TAP_FACTRL ( N , STATUS )

Purpose:
   Return the value of factorial N,  i.e. N! as a floating point number

Language:
   Starlink Fortran 77

Invocation:
   DOUBLE PRECISION TAP_FACTRL
   RESULT = TAP_FACTRL( N, STATUS )

Package:
   TAP: Theoretical Astrophysics Package

Module type:
   double precision function

Description:
   The factorial function is evaluated explicitly the first time it
   is required for a given N > NTOP. The results of n! for all n <= N 
   are stored in an array, so that repeat calculations are not performed. 
   NTOP, initially 0, is reset to N. If N < NTOP, the value of N! is
   recovered from a look up table.

   The program is only valid for N! < the machine largest number. If
   N! = nmax! > this value, an error value is flagged. TAP_FACTRL is not
   defined for -1 < n < nmax.

Accuracy:
   This function is accurate to at least 14 significant figures
   for N <= 169. It has been tested by comparison with the MAPLE
   symbolic computation package (Char, B.W., Geddes, K.O., 
   Gonnet, G.H., Monagan, M.B., and Watt, S.M., 1988, Watcom 
   Publications, Waterloo, Ontario, Canada).

Arguments:
   N = INTEGER (Given)
      The integer N for which N! is required.
   STATUS = INTEGER (Given and Returned)
      The global status.

Returned Value:
   TAP_FACTRL = DOUBLE PRECISION
      N! as a floating point number

---

DOUBLE PRECISION FUNCTION TAP1_POWERG ( Z, STATUS )

Purpose:
   Evaluate power series expansion for ln Gamma (z)

Invocation:
   Should not be invoked except by TAP_GAMMLN
 
Description:
   The Gamma function is evaluated by means of a power series expansion
   valid for |z|<2

      ln (Gamma(1+z)) = -ln(1+z) + z(1-gamma) + sum(n=2,...)(Cn.z^n)

      Cn = (-1)**n * (Zeta(n)-1) / n

      Zeta(n) is the Riemann Zeta function. 

   The constants Cn have been evaluated to 34 d.p. using the symbolic 
   mathematical computation package MAPLEV...

   It is designed to converge to machine precision, but this depends on
   the number of terms (Cn) available. In the current version, the
   TAP parameter TAP__DSMAL is used as convergence criterion. In the
   range 1 < z < 2, the function is accurate to at least 10 decimal 
   places in a comparison with the tables in Abramowitz and Stegun (1971). 

   In the range 0.5 < z < 1.5, the series converges to 1 part in 10^10
   after 15 terms. For z = 2, approximately 30 terms are required.

Arguments:
   Z = DOUBLE PRECISION (Given)
      The real number Z for which ln(GAMMA(X)) is required.
   STATUS = INTEGER (Given and returned)
      The inheritied global status.

Returned Value:
   TAP1_POWERG = DOUBLE PRECISION
      LOG (GAMMA(Z)) as a real number

Notes:
   This function subprogram must _only_ be used in association with
   TAP_GAMMLN or an equivalent front end, in order to check for
   argument validity.

---