HELP GAMMA David Young, July 1987
Computation of the Gamma function for real and complex arguments.
gammaln(<num1>) -> <num2>
gamma(<num1>) -> <num2>
gammainv(<num1>) -> <num2>
Keywords: gamma function, transcendental functions, complex numbers,
numerical methods, factorials.
CONTENTS - (Use <ENTER> g to access required sections)
-- Introduction
-- Examples
-- Accuracy
-- Previous implementation
-- References
-- Related Documentation-- Introduction -------------------------------------------------------
The gamma function is defined as
infinity (z-1) -t
gamma(z) = integral { t * e } dt
0
For positive integer arguments, gamma(n+1) is equal to n!, where n! is
the factorial function of n, that is n*(n-1)*(n-2)*...*3*2.
For real arguments, the function can be found in standard mathematical
libraries. However, for complex arguments an implementation can be hard
to come by. LIB *GAMMA provides an approximation to the function for all
values of the argument. It is also useful for approximating the
factorial function when absolute accuracy is not essential and the
result is within the range of representation of real numbers. The
algorithm is adapted from Press et al. (1988) (see references below).
Procedure GAMMA returns an approximation to the gamma function defined
above.
Procedure GAMMALN returns an approximation to the logarithm of the gamma
function, and should be used whenever a ratio of gamma functions is
needed, to avoid overflow (e.g. in working out binomial coefficients).
Procedure GAMMAINV returns an approximation to the reciprocal of the
function, ie 1/gamma(z). It should be used directly whenever the
reciprocal is required, in order to avoid possible overflow near the
poles of the function.
-- Examples -----------------------------------------------------------
1. Find the factorial of 20.
true -> popdprecision;
gamma(20 + 1) =>
** 2432902007680319000.0
The exact value is in fact 2432902008176640000.
2. Find the logarithm of the gamma function of 2.3 + 2.7i and check it
against Abramowitz & Stegun (1964):
true -> popdprecision;
gammaln(2.3 +: 2.7) =>
** -1.422992_+:2.257905
The table gives the logarithm of the gamma function of 1.3 + 2.7i as
-2.52049157 + 1.13583190i. Using the logarithm of the recurrence formula
z*gamma(z) = gamma(z+1) allows a check on the result:
(-2.52049157 +: 1.13583190) + log(1.3 +: 2.7) =>
** -1.422992_+:2.257905
(If you try this for large values of the imaginary part, bear in mind
that you may need to add multiples of 2*pi to the imaginary part of the
complex logarithm to get agreement.)
-- Accuracy -----------------------------------------------------------
It is recommended that the procedures be run with *POPDPRECISION set
TRUE.
According to Press et al. (1988), the relative accuracy of the results
should be better than about 2 parts in 10^-10 if the real part of the
argument is greater than 1. Other values are calculated using the
reflection formula, and should generally achieve this accuracy except
very close to the poles at 0 and the real negative integers, where the
logarithm of the sine of the argument has to be taken.
The procedures produce results in good agreement with tabulated values
of Abramowitz & Stegun and also agree to better than 1 part in 10^-5
with results obtained using a completely different algorithm based on
formulae given by the same Abramowitz & Stegun. (This is based on a
large number of numerical tests using random points generated over the
complex plane for |z| < 10.) If values for real arguments only are
required, routines from standard libraries such as the the NAG library
may possibly be more reliable and more efficient. If accuracy is
critical the user is advised to carry out an independent test of
accuracy in the context of his or her application.
The accuracy for large positive integers can be tested easily by
comparison with exact calculations of the factorial function.
-- Previous implementation --------------------------------------------
Previous versions of POPLOG included a library based on methods from
Abramowitz & Stegun. This has been replaced the method of Press et al.
(actually derived by Lanczos) because it involves fewer numerical
coefficients, has more accessible documentation, and involves fewer
branches and hence fewer possible discontinuities. It is also a good
idea to find the logarithm, given the large values that arise for many
arguments.
The procedures GAMMA1 and GAMMAINV0 are withdrawn.
-- References ---------------------------------------------------------
Abramowitz, M. & Stegun, I.A. (1964). Handbook of Mathematical Functions
with Formulas, Graphs and Mathematical Tables (Seventh printing, May
1968, with corrections). Washington: National Bureau of Standards.
Press, W.H., Flannery, B.P., Teukolsky, S.A. & Vetterling, W.T. (1988).
Numerical Recipes in C. Cambridge: CUP.
-- Related Documentation ----------------------------------------------
REF *NUMBERS
--- C.all/help/gamma ---------------------------------------------------
--- Copyright University of Sussex 1987. All rights reserved. ----------