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.  The functions described in this section are declared in the header file gsl_sf_gamma.h.

 -- 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 can be reconstructed using the relation \Gamma(x) = sgn * \exp(resultlg).

 -- 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), subject to a and a+x not being negative integers or zero. The Pochhammer symbol is also known as the Apell symbol and sometimes written as (a,x).

 -- 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)) for a > 0, a+x > 0.

 -- 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|) and sgn = \sgn((a)_x) where (a)_x = \Gamma(a + x)/\Gamma(a), subject to a, a+x not being negative integers or zero.

 -- 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).