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

.. code-block:: C++

   #include <gnstlib.hpp>

==============
Main functions
==============

GNSTLIB::gamma
--------------

.. cpp:function:: double gamma(const double x, int& err_id)

   Computes the Gamma function :math:`\Gamma(x)` for real :math:`x`.
   The Gamma function is defined by Euler's integral

   .. math::
      \Gamma(z) = \int_0^{\infty} t^{z - 1} e^{-t} \mathop{dt}, \quad \Re(z) > 0,
   
   and when :math:`\Re(z) < 0`, :math:`\Gamma(z)` is defined by analytic continuation. Important relations are

   .. math:: 
      \Gamma(z + 1) = z \Gamma(z), \quad \text{(recursion)},

   and 

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

   #. **Arguments:**

      * **x**: (double) - *argument*.
         *Constraint*: **x** :math:`\notin \mathbb{Z}^-_0`.

      * **err_id**: (int) - *error identifier*.

   #. **Errors and Warnings:**

      * ``err_id`` = 1: 
         Argument **x** is either too large or too close to 0.0. 
         More precisely **x** :math:`\gtrsim` 171.62437... or **x** :math:`\in` (0.0, 1.0/``MAXDOUBLE``). 
         Result is too large to be represented and :math:`\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, :math:`\infty` is returned. For **x** :math:`= -\infty`, NaN is returned.
      * ``err_id`` = -1: 
         Argument **x** is NaN or an unexpected error occurred. Contact the authors.

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

   Computes the Gamma function :math:`\Gamma(z)` for complex :math:`z`.

   #. **Arguments:**

      * **z**: (complex) - *argument*.
         *Constraint*: **z** :math:`\notin \mathbb{Z}^-_0`.

      * **err_id**: (int) - *error identifier*.

   #. **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 :math:`\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, :math:`\infty` is returned. For **z** :math:`= -\infty`, NaN is returned.
      * ``err_id`` = -1: 
         Argument **z** is NaN or an unexpected error occurred. Contact the authors.

GNSTLIB::gammaln
----------------

.. cpp:function:: double gammaln(const double x, int& err_id)

   Computes the natural logarithm of the absolute value of the Gamma function :math:`\Gamma(x)` for real :math:`x`. For :math:`x < 0`, this function returns the real part of :math:`\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 :cpp:func:`gamma()` and :cpp:func:`gammasign()`.

   #. **Arguments:**

      * **x**: (double) - *argument*.
         *Constraint*: **x** :math:`\notin \mathbb{Z}^-_0`.

      * **err_id**: (int) - *error identifier*.


   #. **Errors and Warnings:**

      * ``err_id`` = 1: 
         Argument **x** is either too large or too close to 0.0. 
         More precisely **x** :math:`\gtrsim` 2.559983...e305 or **x** :math:`\in` (0.0, 1.0/``MAXDOUBLE``). 
         Result is too large to be represented and :math:`\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, :math:`\infty` is returned. For **x** :math:`= -\infty`, :math:`\infty` is returned.
      * ``err_id`` = -1: 
         Argument **x** is NaN or an unexpected error occurred. Contact the authors.

GNSTLIB::gammasign
------------------

.. cpp:function:: int gammasign(const double x)

   Returns the sign of the Gamma function.

   #. **Arguments:**

      * **x**: (double) - *argument*.


GNSTLIB::loggamma
-----------------

.. cpp:function:: 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 :math:`z`. An asymptotic expansion of the logarithm of the gamma function (Stirling' series, see :cpp:func:`stirling()`) is given by

   .. math::
      \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 :math:`R_N(z)` is the error term. The :math:`\log \Gamma(z)` is analytic except for its single branch cut on the negative real axis. The implementation is based on [1].

   #. **Arguments:**

      * **z**: (complex) - *argument*.
         *Constraint*: **z** :math:`\notin \mathbb{Z}^-_0`.

      * **err_id**: (int) - *error identifier*.   

   #. **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 :math:`\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, :math:`\infty` is returned. For **z** :math:`= -\infty`, :math:`\infty` is returned.
      * ``err_id`` = -1: 
         Argument **z** is NaN or an unexpected error occurred. Contact the authors.

GNSTLIB::rgamma
---------------

.. cpp:function:: std::complex<double> rgamma(const std::complex<double> z)

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

   #. **Arguments:**

      * **z**: (complex) - *argument*.

GNSTLIB::qgamma
---------------

.. cpp:function:: double qgamma(const double x, const double y, int& err_id)

   Computes the quotients of two gamma functions, :math:`\frac{\Gamma(a)}{\Gamma(b)}` for :math:`(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.

   #. **Arguments:**

      * **x**: (double) - *argument*.
      * **y**: (double) - *argument*.
      * **err_id**: (int) - *error identifier*.

   #. **Errors and Warnings:**

      * ``err_id`` = 1: 
         Result is too large to be represented and :math:`\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
---------------

.. cpp:function:: double auxgam(const double x, int& err_id)

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

   #. **Arguments:**

      * **x**: (double) - *argument*.
         *Constraint*: **x** :math:`\in [-1,1]`.

      * **err_id**: (int) - *error identifier*.  

   #. **Errors and Warnings:**

      * ``err_id`` = 1: 
         Result is too large to be represented and :math:`\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
------------------

.. cpp:function:: double gammastar(const double x, int& err_id)

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

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

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

   #. **Arguments:**

      * **x**: (double) - *argument*.
         *Constraint*: **x** :math:`>` 0.0.

      * **err_id**: (int) - *error identifier*.

   #. **Errors and Warnings:**

      * ``err_id`` = 1: 
         Argument **x** is too close to 0.0.
         Result is too large to be represented and :math:`\infty` is returned.
      * ``err_id`` = 4:
         Argument **x** :math:`\le` 0.0, NaN is returned.
      * ``err_id`` = -1: 
         Argument **x** is NaN or an unexpected error occurred. Contact the authors.

GNSTLIB::stirling
-----------------

.. cpp:function:: double stirling(const double x, int& err_id)

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

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

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


   #. **Arguments:**

      * **x**: (double) - *argument*.
         *Constraint*: **x** :math:`>` 0.0.

      * **err_id**: (int) - *error identifier*.

   #. **Errors and Warnings:**

      * ``err_id`` = 1: 
         Argument **x** is too close to 0.0. More precisely **x** :math:`\in` (0.0, 1.0/``MAXDOUBLE``).
         Result is too large to be represented and :math:`\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** :math:`\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
------------------

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

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

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

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


GNSTLIB::gammaln_vec
--------------------

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

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


GNSTLIB::loggamma_vec
---------------------

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

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


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

.. cpp:function:: 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.