****************************************************
Exponential, Logarithmic and Trigonometric Integrals
****************************************************

GNSTLIB implements the exponential, logarithmic and trigonometric integrals for real and complex arguments. For inverse functions only
real arguments are supported. This set of functions is included in the main header file:

.. code-block:: C++

   #include <gnstlib.hpp>

=====================
Exponential integrals
=====================

GNSTLIB::ei
-----------
.. cpp:function:: double ei(const double x, int& err_id)

   Computes the exponential integral :math:`\text{Ei}(x)` for real :math:`x`.
   The exponential integral is defined by

   .. math::
      \text{Ei}(z) = \int_{-\infty}^z \frac{e^{-t}}{t} \mathop{dt}, \quad |\text{arg}(-z)| < \pi.

   Specific methods are considered in the vicinity of the simple zero of :math:`\text{Ei}(x_0)`, :math:`x_0 \approx 0.372507`, in order to maintain relative accuracy.

   #. **Arguments:**

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

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

   #. **Errors and Warnings:**

      * ``err_id`` = 1: 
         Argument **x** is too large, more precisely **x** :math:`\gtrsim` 716.351. Result is too large to be represented and :math:`\infty` is returned.
      * ``err_id`` = 2: 
         Argument **x** is too large and negative. 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:`=` 0.0 and :math:`-\infty` is returned.
      * ``err_id`` = -1:
         Argument **x** is NaN or an unexpected error occurred. Contact the authors.

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

   Computes the exponential integral :math:`\text{Ei}(z)` for complex :math:`z`.

   #. **Arguments:**

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

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

   #. **Errors and Warnings:**

      * ``err_id`` = 1: 
         Argument **z** is too large. Result is too large to be represented and :math:`\infty` is returned.
      * ``err_id`` = 2: 
         Argument **z** is too large and negative. 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** :math:`=` 0.0 and :math:`-\infty` is returned.
      * ``err_id`` = -1:
         Argument **z** is NaN or an unexpected error occurred. Contact the authors.

GNSTLIB::e1
-----------

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

   Computes the exponential integral :math:`\text{E}_1(x)` for real :math:`x`, defined by

   .. math::
      \text{E}_1(z) = -\text{Ei}(-z) = \int_z^{\infty} \frac{e^{-t}}{t} \mathop{dt}, \quad |\text{arg}\; z| < \pi.

   :math:`\text{E}_1(z)` is analytic on the plane with a cut along the negative real axis.

   #. **Arguments:**

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

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

   #. **Errors and Warnings:**

      * ``err_id`` = 1: 
         Argument **x** is too large and negative, more precisely **x** :math:`\lesssim` -716.351. Result is too large to be represented and :math:`-\infty` is returned.
      * ``err_id`` = 2: 
         Argument **x** is too large and positive. 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:`=` 0.0 and :math:`\infty` is returned.
      * ``err_id`` = 100:
         Argument **x** :math:`<` 0.0 and only the real part of the principal value of the integral is returned. Note that :math:`\text{E}_1(-x \pm i0) = -\text{Ei}(x) \mp i\pi`, which is implemented if ``z=std::complex<double>(x, 0)`` (**x** :math:`<` 0.0) is passed as input. See example below.
      * ``err_id`` = -1: 
         Argument **x** is NaN or an unexpected error occurred. Contact the authors.


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

   Computes the exponential integral :math:`\text{E}_1(z)` for complex :math:`z`.

   #. **Arguments:**

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

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

   #. **Errors and Warnings:**

      * ``err_id`` = 1: 
         Argument **z** is too large and negative, more precisely :math:`\Re(z)` :math:`\lesssim` -716.351. Result is too large to be represented and complex :math:`-\infty` is returned.
      * ``err_id`` = 2: 
         Argument **z** is too large and :math:`\Re(z) > 0`. 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** :math:`=` 0.0 and complex :math:`\infty` is returned.
      * ``err_id`` = -1: 
         Argument **z** is NaN or an unexpected error occurred. Contact the authors.

   #. **Example:**

      .. code-block:: C++

         #include <iomanip>
         #include <iostream>

         #include <gnstlib.hpp>

         void test_exp()
         {
            int err_id;

            double x1 = 2.2;
            double x2 = -2.2;
            std::complex<double> z2 = std::complex<double>(-2.2, 0.0);

            double r1 = GNSTLIB::e1(x1, err_id);
            std::cout << "E1(x1) = " << std::setprecision(16) << r1 << std::endl;
            std::cout << "err_id = " << err_id << std::endl;

            double r2 = GNSTLIB::e1(x2, err_id);
            std::cout << "E1(x2) = " << std::setprecision(16) << r1 << std::endl;
            std::cout << "err_id = " << err_id << std::endl;

            std::complex<double> rz = GNSTLIB::e1(z2, err_id);
            std::cout << "E1(z2) = " << std::setprecision(16) << r1 << std::endl;
            std::cout << "err_id = " << err_id << std::endl;
         }

      The above example produces the following output::

         E1(x1) = 0.03719113705193257
         err_id = 0
         E1(x2) = -5.732614699814381
         err_id = 100
         E1(z2) = (-5.732614699814381,-3.141592653589793)
         err_id = 0

====================
Logarithmic integral
====================

GNSTLIB::li
-----------

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

   Computes the logarithmic integral :math:`\text{li}(x)` for :math:`x \ge 0`. :math:`\text{li}(x)` is defined by 

   .. math::
      \text{li}(z) = \int_0^z \frac{1}{\log(z)} \mathop{dt}, 

   and can be represented in terms of the exponential integral :math:`\text{Ei}(x)` using the functional relation :math:`\text{li}(x) = \text{Ei}(\log(x))` for :math:`x \neq 1`, which corresponds to a singularity.

   #. **Arguments:**

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

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

   #. **Errors and Warnings:**

      * ``err_id`` = 1: 
         Argument **x** is too large, more precisely **x** :math:`>` ``MAXDOUBLE``. Result is too large to be represented and :math:`\infty` is returned.
      * ``err_id`` = 2: 
         Argument **x** is too close to 0.0. 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:`<` 0.0, the result is complex and it is returned as NaN.
      * ``err_id`` = 5:
         Argument **x** :math:`=` 1.0 and :math:`-\infty` is returned.
      * ``err_id`` = -1: 
         Argument **x** is NaN or an unexpected error occurred. Contact the authors.

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

   Computes the logarithmic integral :math:`\text{li}(z)` for complex :math:`z`.

   #. **Arguments:**

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

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

   #. **Errors and Warnings:**

      * ``err_id`` = 1: 
         Argument **z** is too large. Result is too large to be represented and :math:`\infty` is returned.
      * ``err_id`` = 2: 
         Argument **z** is too close to 0.0. Result is too small to be represented and 0.0 is returned.
      * ``err_id`` = 3:
         The result is a denormalized / subnormal number.
      * ``err_id`` = 5:
         Argument **z** :math:`=` 1.0 and :math:`-\infty` is returned.
      * ``err_id`` = -1: 
         Argument **z** is NaN or an unexpected error occurred. Contact the authors.

=======================
Trigonometric integrals
=======================

GNSTLIB::ci
-----------

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

   Computes the cosine integral :math:`\text{Ci}(x)` for real :math:`x`. :math:`\text{Ci}(x)` is defined for complex :math:`z` by 

   .. math::
      \begin{align}
      \text{Ci}(z) &= -\int_z^{\infty} \frac{\cos(t)}{t} \mathop{dt}, \quad |\text{arg}(-z)| < \pi, \\
      & = \int_0^z \frac{\cos(t) - 1}{t} \mathop{dt} + \log(z) + \gamma.
      \end{align}

   where the path of integration lies on the plane cut along the negative real axis, not passing through the origin.

   #. **Arguments:**

      * **x**: (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: 
         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:`=` 0.0 and :math:`-\infty` is returned.
      * ``err_id`` = 100:
         Argument **x** :math:`<` 0.0 and only the real part of the principal value of the integral is returned. Note that :math:`\text{Ci}(-x) = \text{Ci}(x) + i\pi`, which is implemented if ``z=std::complex<double>(x, 0)`` (**x** :math:`<` 0.0) is passed as input. See similar case in example in :cpp:func:`e1()`.
      * ``err_id`` = -1: 
         Argument **x** is NaN or an unexpected error occurred. Contact the authors.

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

   Compute the cosine integral :math:`\text{Ci}(z)` for complex :math:`z`. 

   #. **Arguments:**

      * **z**: (complex) - *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: 
         Argument **z** 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 **z** :math:`=` 0.0 and :math:`-\infty` is returned.
      * ``err_id`` = -1: 
         Argument **z** is NaN or an unexpected error occurred. Contact the authors.

GNSTLIB::si
-----------

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

   Computes the sine integral :math:`\text{Si}(x)` for real :math:`x`. :math:`\text{Si}(z)` is an odd entire function defined for complex :math:`z` by

   .. math::
      \text{Si}(z) = \int_0^z \frac{\sin(t)}{t} \mathop{dt}.


   #. **Arguments:**

      * **x**: (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`` = -1: 
         Argument **x** is NaN or an unexpected error occurred. Contact the authors.

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

   Compute the sine integral :math:`\text{Si}(z)` for complex :math:`z`.

   #. **Arguments:**

      * **z**: (complex) - *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`` = -1: 
         Argument **z** is NaN or an unexpected error occurred. Contact the authors.

GNSTLIB::sici
-------------

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

   Computes the sine and cosine integral at once for complex :math:`z`.

   #. **Arguments:**

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

      * **si**: (complex) - *sine integral*.

      * **ci**: (complex) - *cosine integral*.
      
      * **err_id**: (int) - *error identifier*.

====================
Hyperbolic integrals
====================

GNSTLIB::chi
------------

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

   Computes the hyperbolic cosine integral :math:`\text{Chi}(x)` for real :math:`x`. :math:`\text{Chi}(x)` is defined in analogy with the cosine integral for complex :math:`z` by

   .. math::
      \text{Chi}(z) = -\int_z^{\infty} \frac{\cosh(t)}{t} \mathop{dt}, \quad |\text{arg}(-z)| < \pi.

   #. **Arguments:**

      * **x**: (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** :math:`=` 0.0 and :math:`-\infty` is returned.
      * ``err_id`` = 100:
         Argument **x** :math:`<` 0.0 and only the real part of the principal value of the integral is returned. Note that :math:`\text{Chi}(-x) = \text{Chi}(x) + i\pi`, which is implemented if ``z=std::complex<double>(x, 0)`` (**x** :math:`<` 0.0) is passed as input. See similar case in example in :cpp:func:`e1()`.
      * ``err_id`` = -1: 
         Argument **x** is NaN or an unexpected error occurred. Contact the authors.

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

   Compute the hyperbolic cosine integral :math:`\text{Chi}(z)` for complex :math:`z`. 

   #. **Arguments:**

      * **z**: (complex) - *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: 
         Argument **z** 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 **z** :math:`=` 0.0 and :math:`-\infty` is returned.
      * ``err_id`` = -1: 
         Argument **z** is NaN or an unexpected error occurred. Contact the authors.

GNSTLIB::shi
------------

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

   Computes the hyperbolic sine integral :math:`\text{Shi}(x)` for real :math:`x`. :math:`\text{Shi}(z)` is defined in analogy with the sine integral for complex :math:`z` by

   .. math::
      \text{Si}(z) = \int_0^z \frac{\sinh(t)}{t} \mathop{dt}.

   #. **Arguments:**

      * **x**: (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`` = -1: 
         Argument **x** is NaN or an unexpected error occurred. Contact the authors.

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

   Compute the hyperbolic sine integral :math:`\text{Shi}(z)` for complex :math:`z`.

   #. **Arguments:**

      * **z**: (complex) - *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`` = -1: 
         Argument **z** is NaN or an unexpected error occurred. Contact the authors.

=================
Inverse functions
=================

GNSTLIB::invei
--------------

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

   Computes inverse of exponential integral :math:`\text{Ei}^{-1}(x)` for :math:`x \in \mathbb{R}`. Numerical inversion is computed by using Halley's method with suitable starting values, some of them based on [1]. Normally, one or two iterations are required.

   #. **Arguments:**

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

   #. **Errors and Warnings:**

      * ``err_id`` = 1: 
         Argument **x** is too large, more precisely **x** :math:`>` ``MAXDOUBLE``. Result is too large to be represented and :math:`\infty` is returned.
      * ``err_id`` = 2: 
         Argument **x** is too large and negative. Result is too small to be represented and 0.0 is returned.
      * ``err_id`` = 3:
         The result is a denormalized / subnormal number.
      * ``err_id`` = -1: 
         Argument **x** is NaN or an unexpected error occurred. Contact the authors.


GNSTLIB::invli
--------------

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

   Computes inverse of logarithmic integral :math:`\text{li}^{-1}(x)` for :math:`x \in \mathbb{R}`. Numerical inversion is computed by using Halley's method with suitable starting values based on https://math.stackexchange.com/questions/853178/numerical-inverse-of-logarithmic-integral. Normally, one or two iterations are required.

   #. **Arguments:**

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

   #. **Errors and Warnings:**

      * ``err_id`` = 1: 
         Argument **x** is too large, more precisely **x** :math:`>` ``MAXDOUBLE`` / 1000. Result is too large to be represented and :math:`\infty` is returned.
      * ``err_id`` = 3:
         The result is a denormalized / subnormal number.
      * ``err_id`` = -1: 
         Argument **x** is NaN or an unexpected error occurred. Contact the authors.

====================
Vectorized functions
====================

GNSTLIB::ei_vec
---------------

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

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

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

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

GNSTLIB::e1_vec
---------------

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

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

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

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

GNSTLIB::li_vec
---------------

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

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

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

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

GNSTLIB::ci_vec
---------------

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

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

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

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

GNSTLIB::si_vec
---------------

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

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

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

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

GNSTLIB::chi_vec
----------------

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

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

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

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

GNSTLIB::shi_vec
----------------

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

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

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

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

==========
References
==========

.. [1] P. Pecina (1984). On the function inverse to the exponential integral function. Bulletin of the Astronomical Institutes of Czechoslovakia.