Incomplete Gamma and Generalized Exponential Integral

GNSTLIB implements normalized and non-normalized incomplete gamma functions and the generalized exponential integral for real parameters and argument. This set of functions is included in the main header file:

#include <gnstlib.hpp>

Incomplete gamma functions

The incomplete gamma functions are defined by

\[\gamma(a, x) = \int_0^x t^{a-1} e^{-t} \mathop{dt}, \quad \Gamma(a, x) = \int_x^{\infty} t^{a-1} e^{-t} \mathop{dt}.\]

Let us define the normalized forms

\[P(a, x) = \frac{1}{\Gamma(a)} \gamma(a, x), \quad Q(a, x) = \frac{1}{\Gamma(a)} \Gamma(a, x),\]

where we assume that \(a\) and \(x\) are real positive numbers. The functions \(P(a, x)\) and \(Q(a, x)\) are the central gamma distribution function and its complementary function, respectively. The central gamma distribution functions satisfy the complementary relation

\[P(a, x) + Q(a, x) = 1.\]

These distribution functions which appear in many problems of applied probability include, as particular case, the standard chi-square probability functions \(P(\chi^2| \nu)\) and \(Q(\chi^2| \nu)\) with parameters \(a = \nu / 2\) and \(x = \chi^2 / 2\). The implementation of the central gamma distribution functions or normalized incomplete gamma functions is described in [1] and is based on GammaCHI package, see http://www.sciencedirect.com/science/article/pii/S0010465515000107.

GNSTLIB::gammainc_p

double gammainc_p(const double a, const double x, int &err_id)

Computes the normalized lower incomplete gamma function \(P(a, x)\) for positive real parameter and argument.

  1. Arguments:

    • a: (double) - argument.

      Constraint: a \(> 0\).

    • x: (double) - argument.

      Constraint: x \(\ge 0\).

    • 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:

      Parameter a \(\le\) 0.0, NaN is returned.

    • err_id = 5:

      Argument x is negative, NaN is returned.

    • err_id = -1:

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

GNSTLIB::gammainc_q

double gammainc_q(const double a, const double x, int &err_id)

Computes the normalized upper incomplete gamma function \(Q(a, x)\) for positive real parameter and argument.

  1. Arguments:

    • a: (double) - argument.

      Constraint: a \(> 0\).

    • x: (double) - argument.

      Constraint: x \(\ge 0\).

    • 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:

      Parameter a \(\le\) 0.0, NaN is returned.

    • err_id = 5:

      Argument x is negative, NaN is returned.

    • err_id = -1:

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

GNSTLIB::gammainc_lower

double gammainc_lower(const double a, const double x, int &err_id)

GNSTLIB::gammainc_upper

double gammainc_upper(const double a, const double x, int &err_id)

GNSTLIB::gammainc_tricomi

double gammainc_tricomi(const double a, const double x, int &err_id)

Generalized exponential integral

GNSTLIB::expint

double expint(const double v, const double x, int &err_id)
Computes the generalized exponential integral \(E_{\nu}(x)\) for \(v \ge 0\) and \(x \ge 0\). The generalized exponential integral is defined by
\[E_{\nu}(x) = \int_1^{\infty} e^{-xt} t^{-\nu} \mathop{dt}, \quad \nu \in \mathbb{R},\, x > 0.\]

The generalized exponential integral can also be expressed in terms of the lower and upper incomplete gamma functions gammainc_lower() and gammainc_upper(), respectively. The implemented algorithm is described in [2].

  1. Arguments:

    • v: (double) - argument.

      Constraint: v \(\ge 0\).

    • x: (double) - argument.

      Constraint: x \(\ge 0\).

    • 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:

      Parameter v is negative, NaN is returned.

    • err_id = 5:

      Argument x is negative, NaN is returned.

    • err_id = -1:

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

double expint(const int n, const double e, const double x, int &err_id)

Computes the generalized exponential integral with positive real parameter and argument. For improved precision for \(\nu \approx |\nu|\) and small \(x\), this routine accepts \(\nu\) split into an integral and a decimal fractional part. Specifically \(\nu = n + e\), where \(|e| \le 0.5\).

  1. Arguments:

    • n: (int) - argument.

      Constraint: v \(\ge 0\).

    • e: (double) - argument.

      Constraint: \(|\) e \(|\le 0.5\).

    • x: (double) - argument.

      Constraint: x \(\ge 0\).

    • 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:

      Parameter n + e is negative, NaN is returned.

    • err_id = 5:

      Argument x is negative, NaN is returned.

    • err_id = 6:

      Argument \(|\) e \(| > 0.5\), NaN is returned.

    • err_id = -1:

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

References

[1]
  1. Gil, J. Segura, and N. M. Temme. Efficient and accurate algorithm for the computation and inversion of the incomplete gamma function ratios. SIAM J. Scientific Computing, 34(6), (2012).

[2]
  1. Navas-Palencia. Fast and accurate algorithm for the generalized exponential integral \(E_{\nu}(x)\) for positive real order. Numerical Algorithms, (2017).