Gamma Functions

GNSTLIB implements the Gamma function and several related functions. These functions support complex argument unless otherwise specified. Gamma-like functions are included in the main header file:

#include <gnstlib.hpp>

Main functions

GNSTLIB::gamma

double gamma(const double x, int &err_id)

Computes the Gamma function \(\Gamma(x)\) for real \(x\). The Gamma function is defined by Euler’s integral

\[\Gamma(z) = \int_0^{\infty} t^{z - 1} e^{-t} \mathop{dt}, \quad \Re(z) > 0,\]

and when \(\Re(z) < 0\), \(\Gamma(z)\) is defined by analytic continuation. Important relations are

\[\Gamma(z + 1) = z \Gamma(z), \quad \text{(recursion)},\]

and

\[\frac{1}{\Gamma(1 + z)\Gamma(1 - z)} = \frac{\sin \pi z}{\pi z}, \quad \text{(reflection formula)}.\]

1. Arguments:

   • x: (double) - argument.

     Constraint: x \(\notin \mathbb{Z}^-_0\).

   • err_id: (int) - error identifier.

2. Errors and Warnings:

   • err_id = 1:

     Argument x is either too large or too close to 0.0. More precisely x \(\gtrsim\) 171.62437... or x \(\in\) (0.0, 1.0/MAXDOUBLE). Result is too large to be represented and \(\infty\) is returned.

   • err_id = 2:

     Argument x is too large and negative non-integer. Result is too small to be represented and 0.0 is returned.

   • err_id = 3:

     The result is a denormalized / subnormal number.

   • err_id = 4:

     Argument x is a negative integer and corresponds to a Gamma function pole, \(\infty\) is returned. For x \(= -\infty\), NaN is returned.

   • err_id = -1:

     Argument x is NaN or an unexpected error occurred. Contact the authors.

std::complex<double> gamma(const std::complex<double> z, int &err_id)

Computes the Gamma function \(\Gamma(z)\) for complex \(z\).

1. Arguments:

   • z: (complex) - argument.

     Constraint: z \(\notin \mathbb{Z}^-_0\).

   • err_id: (int) - error identifier.

2. Errors and Warnings:

   • err_id = 1:

     Argument z is either too large or too close to 0.0. Result is too large to be represented and \(\infty\) is returned.

   • err_id = 2:

     Argument z is too large and negative non-integer. Result is too small to be represented and 0.0 is returned.

   • err_id = 3:

     The result is a denormalized / subnormal number.

   • err_id = 4:

     Argument z is a negative integer and corresponds to a Gamma function pole, \(\infty\) is returned. For z \(= -\infty\), NaN is returned.

   • err_id = -1:

     Argument z is NaN or an unexpected error occurred. Contact the authors.

GNSTLIB::gammaln

double gammaln(const double x, int &err_id)

Computes the natural logarithm of the absolute value of the Gamma function \(\Gamma(x)\) for real \(x\). For \(x < 0\), this function returns the real part of \(\log(\Gamma(x))\).

This function may be useful to avoid computation with complex numbers when working in logspace on the real axis. In this situation the common relation exp(gammaln(x)) = gammasign(x)*gamma(x) might be applied. See gamma() and gammasign().

1. Arguments:

   • x: (double) - argument.

     Constraint: x \(\notin \mathbb{Z}^-_0\).

   • err_id: (int) - error identifier.

2. Errors and Warnings:

   • err_id = 1:

     Argument x is either too large or too close to 0.0. More precisely x \(\gtrsim\) 2.559983...e305 or x \(\in\) (0.0, 1.0/MAXDOUBLE). Result is too large to be represented and \(\infty\) is returned.

   • err_id = 2:

     Argument x is too large and negative non-integer. Result is too small to be represented and 0.0 is returned.

   • err_id = 3:

     The result is a denormalized / subnormal number.

   • err_id = 4:

     Argument x is a negative integer and corresponds to a Gamma function pole, \(\infty\) is returned. For x \(= -\infty\), \(\infty\) is returned.

   • err_id = -1:

     Argument x is NaN or an unexpected error occurred. Contact the authors.

GNSTLIB::gammasign

int gammasign(const double x)

Returns the sign of the Gamma function.

1. Arguments:

   • x: (double) - argument.

GNSTLIB::loggamma

std::complex<double> loggamma(const std::complex<double> z, int &err_id)

Computes the principal branch of the logarithm of the Gamma function for complex \(z\). An asymptotic expansion of the logarithm of the gamma function (Stirling’ series, see stirling()) is given by

\[\log \Gamma(z) = \bigg( z - \frac{1}{2}\bigg) \log(z) - z + \frac{\log(2\pi)}{2} + \sum_{k = 1}^{N-1} \frac{B_{2k}}{2k (2k - 1) z^{2k-1}} + R_N(z),\]

where \(R_N(z)\) is the error term. The \(\log \Gamma(z)\) is analytic except for its single branch cut on the negative real axis. The implementation is based on [1].

1. Arguments:

   • z: (complex) - argument.

     Constraint: z \(\notin \mathbb{Z}^-_0\).

   • err_id: (int) - error identifier.

2. Errors and Warnings:

   • err_id = 1:

     Argument z is either too large or too close to 0.0. Result is too large to be represented and \(\infty\) is returned.

   • err_id = 2:

     Argument z is too large and negative non-integer. Result is too small to be represented and 0.0 is returned.

   • err_id = 3:

     The result is a denormalized / subnormal number.

   • err_id = 4:

     Argument z is a negative integer and corresponds to a Gamma function pole, \(\infty\) is returned. For z \(= -\infty\), \(\infty\) is returned.

   • err_id = -1:

     Argument z is NaN or an unexpected error occurred. Contact the authors.

GNSTLIB::rgamma

std::complex<double> rgamma(const std::complex<double> z)

Computes the reciprocal gamma function \(\frac{1}{\Gamma(z)}\) for complex \(z\).

1. Arguments:

   • z: (complex) - argument.

GNSTLIB::qgamma

double qgamma(const double x, const double y, int &err_id)

Computes the quotients of two gamma functions, \(\frac{\Gamma(a)}{\Gamma(b)}\) for \((a, b) \in \mathbb{R}\).

The quotient of two gamma functions appears frequently in applications. From a computational point of view, problems can arise when trying to compute directly the ratio of functions in the case when the arguments of both functions are large.

1. Arguments:

   • x: (double) - argument.

   • y: (double) - argument.

   • err_id: (int) - error identifier.

2. Errors and Warnings:

   • err_id = 1:

     Result is too large to be represented and \(\infty\) is returned.

   • err_id = 2:

     Result is too small to be represented and 0.0 is returned.

   • err_id = 3:

     The result is a denormalized / subnormal number.

   • err_id = 4:

     Argument x and/or y is a negative integer and corresponds to a Gamma function pole, NaN is returned.

   • err_id = -1:

     Argument x and/or y is NaN or an unexpected error occurred. Contact the authors.

Auxiliary functions

GNSTLIB::auxgam

double auxgam(const double x, int &err_id)

Computes the auxiliary function \(g(x)\) in \(\frac{1}{\Gamma(1 + x)} = 1 + x (x - 1) g(x)\) for \(x \in [-1, 1]\).

1. Arguments:

   • x: (double) - argument.

     Constraint: x \(\in [-1,1]\).

   • err_id: (int) - error identifier.

2. Errors and Warnings:

   • err_id = 1:

     Result is too large to be represented and \(\infty\) is returned.

   • err_id = 2:

     Result is too small to be represented and 0.0 is returned.

   • err_id = 3:

     The result is a denormalized / subnormal number.

   • err_id = 4:

     Argument x is outside the domain of computation, NaN is returned.

   • err_id = -1:

     Argument x and/or y is NaN or an unexpected error occurred. Contact the authors.

GNSTLIB::gammastar

double gammastar(const double x, int &err_id)

Computes the scaled gamma function or regulated gamma function denoted as \(\Gamma^*(x)\) for \(x > 0\). \(\Gamma^*(x)\) is defined by

\[\Gamma^*(x) = \frac{\Gamma(x)}{\sqrt{2\pi} x^{x-1/2} e^{-x}}.\]

When \(x\) is large, this function can be very important in algorithms where the function \(\Gamma(x)\) is involved because \(\Gamma^*(x) = 1 + O(1/x)\).

1. Arguments:

   • x: (double) - argument.

     Constraint: x \(>\) 0.0.

   • err_id: (int) - error identifier.

2. Errors and Warnings:

   • err_id = 1:

     Argument x is too close to 0.0. Result is too large to be represented and \(\infty\) is returned.

   • err_id = 4:

     Argument x \(\le\) 0.0, NaN is returned.

   • err_id = -1:

     Argument x is NaN or an unexpected error occurred. Contact the authors.

GNSTLIB::stirling

double stirling(const double x, int &err_id)

Computes Stirling series \(S(x)\) for \(x > 0\). \(S(x)\) can be written in terms of an expansion involving Bernoulli numbers

\[S(x) = \sum_{n=0}^{\infty} \frac{B_{2n+2}}{(2n+1)(2n+2)z^{2n+1}}.\]

Inserting the first Bernoulli numbers we arrive at the representation

\[S(x) = \frac{1}{12} - \frac{1}{360x^3} + \frac{1}{1260x^5} - \frac{1}{1680x^7} + O(x^{-9}).\]

1. Arguments:

   • x: (double) - argument.

     Constraint: x \(>\) 0.0.

   • err_id: (int) - error identifier.

2. Errors and Warnings:

   • err_id = 1:

     Argument x is too close to 0.0. More precisely x \(\in\) (0.0, 1.0/MAXDOUBLE). Result is too large to be represented and \(\infty\) is returned.

   • err_id = 2:

     Argument x is too large. Result is too small to be represented and 0.0 is returned.

   • err_id = 3:

     The result is a denormalized / subnormal number.

   • err_id = 4:

     Argument x \(\le\) 0.0, NaN is returned.

   • err_id = -1:

     Argument x is NaN or an unexpected error occurred. Contact the authors.

Vectorized functions

GNSTLIB::gamma_vec

void gamma_vec(std::vector<double> &v, std::vector<double> &r, unsigned short option = 0)

void gamma_vec(const int n, const double *v, double *r, unsigned short option = 0)

void gamma_vec(std::vector<std::complex<double>> &v, std::vector<std::complex<double>> &r, unsigned short option = 0)

void gamma_vec(const int n, const std::complex<double> *v, std::complex<double> *r, unsigned short option = 0)

GNSTLIB::gammaln_vec

void gammaln_vec(std::vector<double> &v, std::vector<double> &r, unsigned short option = 0)

void gammaln_vec(const int n, const double *v, double *r, unsigned short option = 0)

GNSTLIB::loggamma_vec

void loggamma_vec(std::vector<std::complex<double>> &v, std::vector<std::complex<double>> &r, unsigned short option = 0)

void loggamma_vec(const int n, const std::complex<double> *v, std::complex<double> *r, unsigned short option = 0)

GNSTLIB::rgamma_vec

void rgamma_vec(std::vector<std::complex<double>> &v, std::vector<std::complex<double>> &r, unsigned short option = 0)

void rgamma_vec(const int n, const std::complex<double> *v, std::complex<double> *r, unsigned short option = 0)

References

[1]	D.E.G. Hare (1997). Computing the Principal Branch of log-Gamma. Journal of Algorithms, 25(2): 221 – 236.