This following routines compute the gamma and beta functions in their full and incomplete forms, as well as various kinds of factorials. The functions described in this section are declared in the header file gsl_sf_gamma.h.

Gamma Functions
---------------

The Gamma function is defined by the following integral,

     \Gamma(x) = \int_0^\infty dt  t^{x-1} \exp(-t)

It is related to the factorial function by \Gamma(n)=(n-1)! for positive integer n.  Further information on the Gamma function can be found in Abramowitz & Stegun, Chapter 6.

 -- Function: double gsl_sf_gamma (double X)
 -- Function: int gsl_sf_gamma_e (double X, gsl_sf_result * RESULT)
     These routines compute the Gamma function \Gamma(x), subject to x not being a negative integer or zero.  The function is computed using the real Lanczos method. The maximum value of x such that \Gamma(x) is not considered an overflow is given by the macro GSL_SF_GAMMA_XMAX and is 171.0.

 -- Function: double gsl_sf_lngamma (double X)
 -- Function: int gsl_sf_lngamma_e (double X, gsl_sf_result * RESULT)
     These routines compute the logarithm of the Gamma function, \log(\Gamma(x)), subject to x not being a negative integer or zero.  For x<0 the real part of \log(\Gamma(x)) is returned, which is equivalent to \log(|\Gamma(x)|).  The function is computed using the real Lanczos method.

 -- Function: int gsl_sf_lngamma_sgn_e (double X, gsl_sf_result * RESULT_LG, double * SGN)
     This routine computes the sign of the gamma function and the logarithm of its magnitude, subject to x not being a negative integer or zero.  The function is computed using the real Lanczos method.  The value of the gamma function and its error can be reconstructed using the relation \Gamma(x) = sgn * \exp(result\_lg), taking into account the two components of RESULT_LG.

 -- Function: double gsl_sf_gammastar (double X)
 -- Function: int gsl_sf_gammastar_e (double X, gsl_sf_result * RESULT)
     These routines compute the regulated Gamma Function \Gamma^*(x) for x > 0. The regulated gamma function is given by,

          \Gamma^*(x) = \Gamma(x)/(\sqrt{2\pi} x^{(x-1/2)} \exp(-x))
                      = (1 + (1/12x) + ...)  for x \to \infty

     and is a useful suggestion of Temme.

 -- Function: double gsl_sf_gammainv (double X)
 -- Function: int gsl_sf_gammainv_e (double X, gsl_sf_result * RESULT)
     These routines compute the reciprocal of the gamma function, 1/\Gamma(x) using the real Lanczos method.

 -- Function: int gsl_sf_lngamma_complex_e (double ZR, double ZI, gsl_sf_result * LNR, gsl_sf_result * ARG)
     This routine computes \log(\Gamma(z)) for complex z=z_r+i z_i and z not a negative integer or zero, using the complex Lanczos method.  The returned parameters are lnr = \log|\Gamma(z)| and arg = \arg(\Gamma(z)) in (-\pi,\pi].  Note that the phase part (ARG) is not well-determined when |z| is very large, due to inevitable roundoff in restricting to (-\pi,\pi].  This will result in a GSL_ELOSS error when it occurs.  The absolute value part (LNR), however, never suffers from loss of precision.

Factorials
----------

Although factorials can be computed from the Gamma function, using the relation n! = \Gamma(n+1) for non-negative integer n, it is usually more efficient to call the functions in this section, particularly for small values of n, whose factorial values are maintained in hardcoded tables.

 -- Function: double gsl_sf_fact (unsigned int N)
 -- Function: int gsl_sf_fact_e (unsigned int N, gsl_sf_result * RESULT)
     These routines compute the factorial n!.  The factorial is related to the Gamma function by n! = \Gamma(n+1). The maximum value of n such that n! is not considered an overflow is given by the macro GSL_SF_FACT_NMAX and is 170.

 -- Function: double gsl_sf_doublefact (unsigned int N)
 -- Function: int gsl_sf_doublefact_e (unsigned int N, gsl_sf_result * RESULT)
     These routines compute the double factorial n!! = n(n-2)(n-4) \dots. The maximum value of n such that n!! is not considered an overflow is given by the macro GSL_SF_DOUBLEFACT_NMAX and is 297.

 -- Function: double gsl_sf_lnfact (unsigned int N)
 -- Function: int gsl_sf_lnfact_e (unsigned int N, gsl_sf_result * RESULT)
     These routines compute the logarithm of the factorial of N, \log(n!).  The algorithm is faster than computing \ln(\Gamma(n+1)) via gsl_sf_lngamma for n < 170, but defers for larger N.

 -- Function: double gsl_sf_lndoublefact (unsigned int N)
 -- Function: int gsl_sf_lndoublefact_e (unsigned int N, gsl_sf_result * RESULT)
     These routines compute the logarithm of the double factorial of N, \log(n!!).

 -- Function: double gsl_sf_choose (unsigned int N, unsigned int M)
 -- Function: int gsl_sf_choose_e (unsigned int N, unsigned int M, gsl_sf_result * RESULT)
     These routines compute the combinatorial factor n choose m = n!/(m!(n-m)!)

 -- Function: double gsl_sf_lnchoose (unsigned int N, unsigned int M)
 -- Function: int gsl_sf_lnchoose_e (unsigned int N, unsigned int M, gsl_sf_result * RESULT)
     These routines compute the logarithm of n choose m.  This is equivalent to the sum \log(n!) - \log(m!) - \log((n-m)!).

 -- Function: double gsl_sf_taylorcoeff (int N, double X)
 -- Function: int gsl_sf_taylorcoeff_e (int N, double X, gsl_sf_result * RESULT)
     These routines compute the Taylor coefficient x^n / n! for $x \ge 0$ x >= 0, $n \ge 0$ n >= 0.

Pochhammer Symbol
-----------------

 -- Function: double gsl_sf_poch (double A, double X)
 -- Function: int gsl_sf_poch_e (double A, double X, gsl_sf_result * RESULT)
     These routines compute the Pochhammer symbol (a)_x = \Gamma(a + x)/\Gamma(a).  The Pochhammer symbol is also known as the Apell symbol and sometimes written as (a,x).  When a and a+x are negative integers or zero, the limiting value of the ratio is returned.

 -- Function: double gsl_sf_lnpoch (double A, double X)
 -- Function: int gsl_sf_lnpoch_e (double A, double X, gsl_sf_result * RESULT)
     These routines compute the logarithm of the Pochhammer symbol, \log((a)_x) = \log(\Gamma(a + x)/\Gamma(a)).

 -- Function: int gsl_sf_lnpoch_sgn_e (double A, double X, gsl_sf_result * RESULT, double * SGN)
     These routines compute the sign of the Pochhammer symbol and the logarithm of its magnitude.  The computed parameters are result = \log(|(a)_x|) with a corresponding error term, and sgn = \sgn((a)_x) where (a)_x = \Gamma(a + x)/\Gamma(a).

 -- Function: double gsl_sf_pochrel (double A, double X)
 -- Function: int gsl_sf_pochrel_e (double A, double X, gsl_sf_result * RESULT)
     These routines compute the relative Pochhammer symbol ((a)_x - 1)/x where (a)_x = \Gamma(a + x)/\Gamma(a).

Incomplete Gamma Functions
--------------------------

 -- Function: double gsl_sf_gamma_inc (double A, double X)
 -- Function: int gsl_sf_gamma_inc_e (double A, double X, gsl_sf_result * RESULT)
     These functions compute the unnormalized incomplete Gamma Function $\Gamma(a,x) = \int_x^\infty dt\, t^(a-1) \exp(-t)$ \Gamma(a,x) = \int_x^\infty dt t^{a-1} \exp(-t) for a real and $x \ge 0$ x >= 0.

 -- Function: double gsl_sf_gamma_inc_Q (double A, double X)
 -- Function: int gsl_sf_gamma_inc_Q_e (double A, double X, gsl_sf_result * RESULT)
     These routines compute the normalized incomplete Gamma Function $Q(a,x) = 1/\Gamma(a) \int_x^\infty dt\, t^(a-1) \exp(-t)$ Q(a,x) = 1/\Gamma(a) \int_x^\infty dt t^{a-1} \exp(-t) for a > 0, $x \ge 0$ x >= 0.

 -- Function: double gsl_sf_gamma_inc_P (double A, double X)
 -- Function: int gsl_sf_gamma_inc_P_e (double A, double X, gsl_sf_result * RESULT)
     These routines compute the complementary normalized incomplete Gamma Function $P(a,x) = 1 - Q(a,x) = 1/\Gamma(a) \int_0^x dt\, t^(a-1) \exp(-t)$ P(a,x) = 1 - Q(a,x) = 1/\Gamma(a) \int_0^x dt t^{a-1} \exp(-t) for a > 0, $x \ge 0$ x >= 0.

     Note that Abramowitz & Stegun call P(a,x) the incomplete gamma function (section 6.5).

Beta Functions
--------------

 -- Function: double gsl_sf_beta (double A, double B)
 -- Function: int gsl_sf_beta_e (double A, double B, gsl_sf_result * RESULT)
     These routines compute the Beta Function, B(a,b) = \Gamma(a)\Gamma(b)/\Gamma(a+b) subject to a and b not being negative integers.

 -- Function: double gsl_sf_lnbeta (double A, double B)
 -- Function: int gsl_sf_lnbeta_e (double A, double B, gsl_sf_result * RESULT)
     These routines compute the logarithm of the Beta Function, \log(B(a,b)) subject to a and b not being negative integers.

Incomplete Beta Function
------------------------

 -- Function: double gsl_sf_beta_inc (double A, double B, double X)
 -- Function: int gsl_sf_beta_inc_e (double A, double B, double X, gsl_sf_result * RESULT)
     These routines compute the normalized incomplete Beta function I_x(a,b)=B_x(a,b)/B(a,b) where $B_x(a,b) = \int_0^x t^a-1 (1-t)^b-1 dt$ B_x(a,b) = \int_0^x t^{a-1} (1-t)^{b-1} dt for $0 \le x \le 1$ 0 <= x <= 1.   For a > 0, b > 0 the value is computed using a continued fraction expansion.  For all other values it is computed using the relation $I_x(a,b,x) = (1/a) x^a _2F_1(a,1-b,a+1,x)/B(a,b)$ I_x(a,b,x) = (1/a) x^a 2F1(a,1-b,a+1,x)/B(a,b).