The 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 :file:`gsl_sf_gamma.h`.

Gamma Functions

---------------

.. index:: gamma functions

The Gamma function is defined by the following integral,

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

It is related to the factorial function by :math:`\Gamma(n) = (n-1)!`

for positive integer :math:`n`.  Further information on the Gamma function

can be found in Abramowitz & Stegun, Chapter 6.  

.. function:: double gsl_sf_gamma (double x)

              int gsl_sf_gamma_e (double x, gsl_sf_result * result)

   These routines compute the Gamma function :math:`\Gamma(x)`, subject to :math:`x`

   not being a negative integer or zero.  The function is computed using the real

   Lanczos method. The maximum value of :math:`x` such that :math:`\Gamma(x)` is not

   considered an overflow is given by the macro :macro:`GSL_SF_GAMMA_XMAX`

   and is 171.0.

.. exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EROUND

.. index:: logarithm of Gamma function

.. function:: double gsl_sf_lngamma (double x)

              int gsl_sf_lngamma_e (double x, gsl_sf_result * result)

   These routines compute the logarithm of the Gamma function,

   :math:`\log(\Gamma(x))`, subject to :math:`x` not being a negative

   integer or zero.  For :math:`x < 0` the real part of :math:`\log(\Gamma(x))` is

   returned, which is equivalent to :math:`\log(|\Gamma(x)|)`.  The function

   is computed using the real Lanczos method.

.. exceptions: GSL_EDOM, GSL_EROUND

.. 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 :math:`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 

   :math:`\Gamma(x) = sgn * \exp(result\_lg)`, taking into account the two 

   components of :data:`result_lg`.

.. exceptions: GSL_EDOM, GSL_EROUND

.. index:: Regulated Gamma function

.. function:: double gsl_sf_gammastar (double x)

              int gsl_sf_gammastar_e (double x, gsl_sf_result * result)

   These routines compute the regulated Gamma Function :math:`\Gamma^*(x)`

   for :math:`x > 0`. The regulated gamma function is given by,

   .. only:: not texinfo

      .. math::

         \Gamma^*(x) &= \Gamma(x)/(\sqrt{2\pi} x^{(x-1/2)} \exp(-x))\cr

                     &= \left(1 + {1 \over 12x} + ...\right) \quad\hbox{for~} x\to \infty\cr

   .. only:: texinfo

      ::

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

.. exceptions: GSL_EDOM

.. index:: Reciprocal Gamma function

.. function:: double gsl_sf_gammainv (double x)

              int gsl_sf_gammainv_e (double x, gsl_sf_result * result)

   These routines compute the reciprocal of the gamma function,

   :math:`1/\Gamma(x)` using the real Lanczos method.

.. exceptions: GSL_EUNDRFLW, GSL_EROUND

.. index:: Complex Gamma function

.. function:: int gsl_sf_lngamma_complex_e (double zr, double zi, gsl_sf_result * lnr, gsl_sf_result * arg)

   This routine computes :math:`\log(\Gamma(z))` for complex :math:`z = z_r + i z_i`

   and :math:`z` not a negative integer or zero, using the complex Lanczos

   method.  The returned parameters are :math:`lnr = \log|\Gamma(z)|` and

   :math:`arg = \arg(\Gamma(z))` in :math:`(-\pi,\pi]`.  Note that the phase

   part (:data:`arg`) is not well-determined when :math:`|z|` is very large,

   due to inevitable roundoff in restricting to :math:`(-\pi,\pi]`.  This

   will result in a :macro:`GSL_ELOSS` error when it occurs.  The absolute

   value part (:data:`lnr`), however, never suffers from loss of precision.

.. exceptions: GSL_EDOM, GSL_ELOSS

Factorials

----------

.. index:: factorial

Although factorials can be computed from the Gamma function, using

the relation :math:`n! = \Gamma(n+1)` for non-negative integer :math:`n`,

it is usually more efficient to call the functions in this section,

particularly for small values of :math:`n`, whose factorial values are

maintained in hardcoded tables.

.. index:: factorial

.. function:: double gsl_sf_fact (unsigned int n)

              int gsl_sf_fact_e (unsigned int n, gsl_sf_result * result)

   These routines compute the factorial :math:`n!`.  The factorial is

   related to the Gamma function by :math:`n! = \Gamma(n+1)`.

   The maximum value of :math:`n` such that :math:`n!` is not

   considered an overflow is given by the macro :macro:`GSL_SF_FACT_NMAX`

   and is 170.

.. exceptions: GSL_EDOM, GSL_EOVRFLW

.. index:: double factorial

.. function:: double gsl_sf_doublefact (unsigned int n)

              int gsl_sf_doublefact_e (unsigned int n, gsl_sf_result * result)

   These routines compute the double factorial :math:`n!! = n(n-2)(n-4) \dots`. 

   The maximum value of :math:`n` such that :math:`n!!` is not

   considered an overflow is given by the macro :macro:`GSL_SF_DOUBLEFACT_NMAX`

   and is 297.

.. exceptions: GSL_EDOM, GSL_EOVRFLW

.. index:: logarithm of factorial

.. function:: double gsl_sf_lnfact (unsigned int n)

              int gsl_sf_lnfact_e (unsigned int n, gsl_sf_result * result)

   These routines compute the logarithm of the factorial of :data:`n`,

   :math:`\log(n!)`.  The algorithm is faster than computing

   :math:`\ln(\Gamma(n+1))` via :func:`gsl_sf_lngamma` for :math:`n < 170`,

   but defers for larger :data:`n`.

.. exceptions: none

.. index:: logarithm of double factorial

.. function:: double gsl_sf_lndoublefact (unsigned int n)

              int gsl_sf_lndoublefact_e (unsigned int n, gsl_sf_result * result)

   These routines compute the logarithm of the double factorial of :data:`n`,

   :math:`\log(n!!)`.

.. exceptions: none

.. index:: combinatorial factor C(m,n)

.. function:: double gsl_sf_choose (unsigned int n, unsigned int m)

              int gsl_sf_choose_e (unsigned int n, unsigned int m, gsl_sf_result * result)

   These routines compute the combinatorial factor :code:`n choose m`

   :math:`= n!/(m!(n-m)!)`

.. exceptions: GSL_EDOM, GSL_EOVRFLW

.. index:: logarithm of combinatorial factor C(m,n)

.. function:: double gsl_sf_lnchoose (unsigned int n, unsigned int m)

              int gsl_sf_lnchoose_e (unsigned int n, unsigned int m, gsl_sf_result * result)

   These routines compute the logarithm of :code:`n choose m`.  This is

   equivalent to the sum :math:`\log(n!) - \log(m!) - \log((n-m)!)`.

.. exceptions: GSL_EDOM 

.. index::

   single: Taylor coefficients, computation of

.. function:: double gsl_sf_taylorcoeff (int n, double x)

              int gsl_sf_taylorcoeff_e (int n, double x, gsl_sf_result * result)

   These routines compute the Taylor coefficient :math:`x^n / n!` for 

   :math:`x \ge 0`, :math:`n \ge 0`

.. exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW

.. _pochhammer-symbol:

Pochhammer Symbol

-----------------

.. index::

   single: Pochhammer symbol

   single:  Apell symbol, see Pochhammer symbol

.. function:: double gsl_sf_poch (double a, double x)

              int gsl_sf_poch_e (double a, double x, gsl_sf_result * result)

   These routines compute the Pochhammer symbol :math:`(a)_x = \Gamma(a + x)/\Gamma(a)`.

   The Pochhammer symbol is also known as the Apell symbol and

   sometimes written as :math:`(a,x)`.  When :math:`a` and :math:`a + x` 

   are negative integers or zero, the limiting value of the ratio is returned. 

.. exceptions:  GSL_EDOM, GSL_EOVRFLW

.. index:: logarithm of Pochhammer symbol

.. function:: double gsl_sf_lnpoch (double a, double x)

              int gsl_sf_lnpoch_e (double a, double x, gsl_sf_result * result)

   These routines compute the logarithm of the Pochhammer symbol,

   :math:`\log((a)_x) = \log(\Gamma(a + x)/\Gamma(a))`.

.. exceptions:  GSL_EDOM

.. 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 :math:`result = \log(|(a)_x|)`

   with a corresponding error term,  and :math:`sgn = \sgn((a)_x)` where :math:`(a)_x = \Gamma(a + x)/\Gamma(a)`.

.. exceptions:  GSL_EDOM

.. index:: relative Pochhammer symbol

.. function:: double gsl_sf_pochrel (double a, double x)

              int gsl_sf_pochrel_e (double a, double x, gsl_sf_result * result)

   These routines compute the relative Pochhammer symbol :math:`((a)_x - 1)/x`

   where :math:`(a)_x = \Gamma(a + x)/\Gamma(a)`.

.. exceptions:  GSL_EDOM

Incomplete Gamma Functions

--------------------------

.. index::

   single: non-normalized incomplete Gamma function

   single: unnormalized incomplete Gamma function

.. function:: double gsl_sf_gamma_inc (double a, double x)

              int gsl_sf_gamma_inc_e (double a, double x, gsl_sf_result * result)

   These functions compute the unnormalized incomplete Gamma Function

   :math:`\Gamma(a,x) = \int_x^\infty dt t^{(a-1)} \exp(-t)`

   for :math:`a` real and :math:`x \ge 0`.

.. exceptions: GSL_EDOM

.. index:: incomplete Gamma function

.. function:: double gsl_sf_gamma_inc_Q (double a, double x)

              int gsl_sf_gamma_inc_Q_e (double a, double x, gsl_sf_result * result)

   These routines compute the normalized incomplete Gamma Function

   :math:`Q(a,x) = 1/\Gamma(a) \int_x^\infty dt t^{(a-1)} \exp(-t)`

   for :math:`a > 0`, :math:`x \ge 0`.

.. exceptions: GSL_EDOM

.. index:: complementary incomplete Gamma function

.. function:: double gsl_sf_gamma_inc_P (double a, double x)

              int gsl_sf_gamma_inc_P_e (double a, double x, gsl_sf_result * result)

   These routines compute the complementary normalized incomplete Gamma Function

   :math:`P(a,x) = 1 - Q(a,x) = 1/\Gamma(a) \int_0^x dt t^{(a-1)} \exp(-t)`

   for :math:`a > 0`, :math:`x \ge 0`.

   Note that Abramowitz & Stegun call :math:`P(a,x)` the incomplete gamma

   function (section 6.5).

.. exceptions: GSL_EDOM

Beta Functions

--------------

.. index:: Beta function

.. function:: double gsl_sf_beta (double a, double b)

              int gsl_sf_beta_e (double a, double b, gsl_sf_result * result)

   These routines compute the Beta Function, :math:`B(a,b) = \Gamma(a)\Gamma(b)/\Gamma(a+b)`

   subject to :math:`a` and :math:`b` not being negative integers.

.. exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW

.. index:: logarithm of Beta function

.. function:: double gsl_sf_lnbeta (double a, double b)

              int gsl_sf_lnbeta_e (double a, double b, gsl_sf_result * result)

   These routines compute the logarithm of the Beta Function, :math:`\log(B(a,b))`

   subject to :math:`a` and :math:`b` not being negative integers.

.. exceptions: GSL_EDOM

Incomplete Beta Function

------------------------

.. index::

   single: incomplete Beta function, normalized

   single: normalized incomplete Beta function

   single: Beta function, incomplete normalized 

.. function:: double gsl_sf_beta_inc (double a, double b, double x)

              int gsl_sf_beta_inc_e (double a, double b, double x, gsl_sf_result * result)

   These routines compute the normalized incomplete Beta function

   :math:`I_x(a,b) = B_x(a,b) / B(a,b)` where

   .. math:: B_x(a,b) = \int_0^x t^{a-1} (1-t)^{b-1} dt

   for :math:`0 \le x \le 1`.

   For :math:`a > 0`, :math:`b > 0` the value is computed using

   a continued fraction expansion.  For all other values it is computed using 

   the relation

   .. only:: not texinfo

      .. math:: I_x(a,b,x) = (1/a) x^a {}_2F_1(a,1-b,a+1,x)/B(a,b)

   .. only:: texinfo

      ::

         I_x(a,b,x) = (1/a) x^a 2F1(a,1-b,a+1,x) / B(a,b)