\input{share/preamble}

  %\VignetteIndexEntry{expint user manual}
  %\VignettePackage{expint}
  %\VignetteDepends{pracma,gsl}
  %\SweaveUTF8

  \title{\pkg{expint}: Exponential integral and incomplete gamma function}
  \author{Vincent Goulet \\ Université Laval}
  \date{}

  %% Additional commands specific to this document
  \newcommand{\Ei}{\operatorname{Ei}}

<<echo=FALSE>>=
library(expint)
options(width = 60)
@

\begin{document}

\maketitle

\section{Introduction}
\label{sec:introduction}

The exponential integral
\begin{equation*}
  E_1(x) = \int_x^\infty \frac{e^{-t}}{t}\, dt,
  \quad x \in \mathbb{R}
\end{equation*}
and the incomplete gamma function
\begin{equation*}
  \Gamma(a, x) = \int_x^\infty t^{a-1} e^{-t}\, dt,
  \quad x > 0, \quad a \in \mathbb{R}
\end{equation*}
are two closely related functions that arise in various
fields of mathematics.

\pkg{expint} is a small package provides facilities to compute the
exponential integral and the incomplete gamma function. Furthermore,
and perhaps most conveniently for R package developers, the package
also gives easy access to the underlying C workhorses through an API.
The C routines are derived from the GNU Scientific Library
\citep[GSL;][]{GSL}.

The package \pkg{expint} started its life in version~2.0-0 of
\pkg{actuar} \citep{actuar}, where I extended the range of admissible
values in the computation of limited expected value functions. This
required an incomplete gamma function that accepts negative values of
argument $a$, as explained at the beginning of Appendix~A of
\citet{LossModels4e}.


\section{Exponential integral}
\label{sec:expint}

\citet[Section~5.1]{Abramowitz:1972} first define the exponential
integral as
\begin{equation}
  \label{eq:E1}
  E_1(x) = \int_x^\infty \frac{e^{-t}}{t}\, dt.
\end{equation}

An alternative definition (to be understood in terms of the Cauchy
principal value due to the singularity of the integrand at zero) is
\begin{equation*}
  \Ei(x) = - \int_{-x}^\infty \frac{e^{-t}}{t}\, dt
         = \int_{-\infty}^x \frac{e^t}{t}\, dt, \quad x > 0.
\end{equation*}
The above two definitions are related as follows:
\begin{equation}
  \label{eq:Ei_vs_E1}
  E_1(-x) = - \Ei(x), \quad x > 0.
\end{equation}

The exponential integral can also generalized to
\begin{equation*}
  E_n(x) = \int_1^\infty \frac{e^{-xt}}{t^n}\, dt, \quad
  n = 0, 1, 2, \dots, \quad x > 0,
\end{equation*}
where $n$ is then the \emph{order} of the integral. The latter
expression is closely related to the incomplete gamma function
(\autoref{sec:gammainc}) as follows:
\begin{equation}
  \label{eq:En_vs_gammainc}
  E_n(x) = x^{n - 1} \Gamma(1 - n, x).
\end{equation}
One should note that the first argument of function $\Gamma$ is
negative for $n > 1$.

The following recurrence relation holds between exponential integrals
of successive orders:
\begin{equation}
  \label{eq:En:recurrence}
  E_{n+1}(x) = \frac{1}{n} [e^{-x} - x E_n(x)].
\end{equation}
Finally, $E_n(x)$ has the following asymptotic expansion:
\begin{equation}
  \label{eq:En:asymptotic}
  E_n(x) \asymp \frac{e^{-x}}{x}
  \left(
    1 - \frac{n}{x} + \frac{n(n+1)}{x^2} - \frac{n(n+1)(n+2)}{x^3} +
    \dots
  \right).
\end{equation}


\section{Incomplete gamma function}
\label{sec:gammainc}

From a probability theory perspective, the incomplete gamma function
is usually defined as
\begin{equation*}
  P(a, x) = \frac{1}{\Gamma(a)} \int_0^x t^{a - 1} e^{-t}\, dt,
  \quad x > 0, \quad a > 0.
\end{equation*}
Function \code{pgamma} already implements this function in
R (just note the differing order of the arguments).

Now, the definition of the incomplete gamma function of interest for
this package is rather the following
\citep[Section~6.5]{Abramowitz:1972}:
\begin{equation}
  \label{eq:gammainc}
  \Gamma(a, x) = \int_x^\infty t^{a-1} e^{-t}\, dt,
  \quad x > 0, \quad a \in \mathbb{R}.
\end{equation}
Note that $a$ can be negative with this definition. Of course, for
$a > 0$ one has
\begin{equation}
  \label{eq:gammainc_vs_pgamma}
  \Gamma(a, x) = \Gamma(a) [1 - P(a, x)].
\end{equation}

Integration by parts of the integral in \eqref{eq:gammainc} yields the
recursive relation
\begin{equation}
  \label{eq:gammainc_recursion}
  \Gamma(a, x) = -\frac{x^a e^{-x}}{a} + \frac{1}{a} \Gamma(a + 1, x).
\end{equation}
When $a < 0$, this relation can be used repeatedly $k$ times until
$a + k$ is a positive number. The right hand side can then be
evaluated with \eqref{eq:gammainc_vs_pgamma}. If
$a = 0, -1, -2, \dots$, this calculation requires the value of
\begin{equation*}
  G(0, x) = \int_x^\infty \frac{e^{-t}}{t}\, dt = E_1(x),
\end{equation*}
the exponential integral defined in \eqref{eq:E1}.


\section{R interfaces}
\label{sec:interfaces}

\pkg{expint} provides one main and four auxiliary R functions to
compute the exponential integral, and one function to compute the
incomplete gamma function. Their signatures are the following:
\begin{Schunk}
\begin{Sinput}
expint(x, order = 1L, scale = FALSE)
expint_E1(x, scale = FALSE)
expint_E2(x, scale = FALSE)
expint_En(x, order, scale = FALSE)
expint_Ei(x, scale = FALSE)
gammainc(a, x)
\end{Sinput}
\end{Schunk}

Let us first go over function \code{gammainc} since there is less to
discuss. The function takes in argument two vectors or real numbers
(non-negative for argument \code{x}) and returns the value of
$\Gamma(a, x)$. The function is vectorized in arguments \code{a} and
\code{x}, so it works similar to, say, \code{pgamma}.

We now turn to the \code{expint} family of functions. The function
\code{expint} is a unified interface to compute exponential integrals
$E_n(x)$ of any (non-negative) order, with default the most common
case $E_1(x)$. The function is vectorized in arguments \code{x} and
\code{order}. In other words, one can compute the exponential integral
of a different order for each value of $x$.
<<echo=TRUE>>=
expint(c(1.275, 10, 12.3), order = 1:3)
@

The argument \code{order} should be a vector of integers. Non-integer
values are silently coerced to integers using truncation towards zero.

When the argument \code{scale} is \code{TRUE}, the result is scaled by
$e^{x}$.

The functions \code{expint\_E1}, \code{expint\_E2} and \code{expint\_En} are
simpler, slightly faster ways to directly compute exponential
integrals $E_1(x)$, $E_2(x)$ and $E_n(x)$, the latter for a \emph{single}
order $n$ (the first value of \code{order} if \code{order} is a
vector).
<<echo=TRUE>>=
expint_E1(1.275)
expint_E2(10)
expint_En(12.3, order = 3L)
@

Finally, the function \code{expint\_Ei} is provided as a convenience to
compute $\Ei(x)$ using \eqref{eq:Ei_vs_E1}.
<<echo=TRUE>>=
expint_Ei(5)
-expint_E1(-5)     # same
@


\section{Accessing the C routines}
\label{sec:api}

The actual workhorses behind the R functions of
\autoref{sec:interfaces} are C routines with the following prototypes:
\begin{Schunk}
\begin{Sinput}
double expint_E1(double x, int scale);
double expint_E2(double x, int scale);
double expint_En(double x, int order, int scale);
double gamma_inc(double a, double x);
\end{Sinput}
\end{Schunk}

\pkg{expint} makes these routines available to other packages through
declarations in the header file \file{include/expintAPI.h} in the
package installation directory. If you want to use a routine --- say
\code{expint\_E1} --- in your package \pkg{pkg}, proceed as follows:
\begin{enumerate}
\item Add the package \pkg{expint} to the \code{Imports} and \code{LinkingTo}
  directives of the \file{DESCRIPTION} file of \pkg{pkg};
\item Add an entry \samp{import(expint)} in the \file{NAMESPACE} file
  of \pkg{pkg};
\item Define the routine with a call to \code{R\_GetCCallable} in the
  initialization routine \code{R\_init\_pkg} of \pkg{pkg}
  \citep[Section~5.4]{WRE}. For the current example, the file
  \file{src/init.c} of \pkg{pkg} would contain the following code:
\begin{Schunk}
\begin{Sinput}
void R_init_pkg(DllInfo *dll)
{
    R_registerRoutines(/* native routine registration */);

    pkg_expint_E1 = (double(*)(double,int,int))
                    R_GetCCallable("expint", "expint_E1");
}
\end{Sinput}
\end{Schunk}
\item Define a native routine interface, say \code{pkg\_expint\_E1} to
  avoid any name clash, in \file{src/init.c} that will call
  \code{expint\_E1}:
\begin{Schunk}
\begin{Sinput}
double(*pkg_expint_E1)(double,int);
\end{Sinput}
\end{Schunk}
\item Declare the routine in a header file of \pkg{pkg} with the
  keyword \code{extern} to expose the interface to all routines of the
  package. In our example, \file{src/pkg.h} would contain:
\begin{Schunk}
\begin{Sinput}
extern double(*pkg_expint_E1)(double,int);
\end{Sinput}
\end{Schunk}
\item Include the package header file \file{pkg.h} in any C file
  making use of the routine \code{pkg\_expint\_E1}.
\end{enumerate}

To help developers get started, \pkg{expint} ships with a complete
test package implementing the above; see the \file{example\_API}
sub-directory in the installation directory. This test package uses
the \code{.External} R to C interface and, as a bonus, shows how to
vectorize an R function on the C side (the code for this being mostly
derived from base R).

There are various ways to define a package API. The approach described
above was derived from the package \pkg{zoo} \citep{zoo}. The package
\pkg{xts} \citep{xts} --- and probably a few others on CRAN --- draws
from \pkg{Matrix} \citep{Matrix} to propose a somewhat simpler
approach where the API exposes routines that can be used directly in a
package. However, the provided header file can be included only once
in a package, otherwise one gets \samp{duplicate symbols} errors at
link time. This constraint does not show in the example provided with
\pkg{xts} or in packages \pkg{RcppXts} \citep{RcppXts} and \pkg{TTR}
\citep{TTR} that link to it (the only two at the time of writing). A
way around the issue is to define a native routine calling the
routines exposed in the API. In this scenario, tests I conducted
proved the approach I retained to be up to 10\% faster most of the
time.


\section{Implementation details}
\label{sec:implementation}

As already stated, the C routines mentioned in
\autoref{sec:api} are derived from code in the GNU Scientific Library
\citep{GSL}.

For exponential integrals, the main routine \code{expint\_E1} computes
$E_1(x)$ using Chebyshev expansions \citep[chapter~3]{Gil:2007}.
Routine \code{expint\_E2} computes $E_2(x)$ using \code{expint\_E1}
with relation \eqref{eq:En:recurrence} for $x < 100$, and using the
asymptotic expression \eqref{eq:En:asymptotic} otherwise. Routine
\code{expint\_En} simply relies on \code{gamma\_inc} to compute
$E_n(x)$ for $n > 2$ through relation \eqref{eq:En_vs_gammainc}.

For the sake of providing routines that better fit within the
R ecosystem and coding style, I made the following changes
to the original GSL code:
\begin{enumerate}
\item routines compute a single value and return their result by
  value;
\item accordingly, calculation of the approximation error is dropped
  in all routines;
\item \code{gamma\_inc} computes $\Gamma(a, x)$ for $a > 0$ with
  \eqref{eq:gammainc_vs_pgamma} using the routines \code{gammafn} and
  \code{pgamma} of the R API, rather than using the GSL routines, as
  the example below illustrates;
<<echo=FALSE>>=
op <- options() # remember default number of digits
@
<<echo=TRUE>>=
options(digits = 20)
gammainc(1.2, 3)
gamma(1.2) * pgamma(3, 1.2, 1, lower = FALSE)
@
<<echo=FALSE>>=
options(op)     # restore defaults
@
\item \code{gamma\_inc} computes $\Gamma(a, x)$ for $-0.5 < a < 0$
  using the recursion \eqref{eq:gammainc_recursion} instead of a
  series expansion as in the GSL routines, thereby relying on the
  accuracy of \code{pgamma} near $a = 0.5$ (fixes
  \href{https://gitlab.com/vigou3/expint/-/issues/2}{issue \#2}; see
  \autoref{sec:computation_gamma_inc} for additional details).
\end{enumerate}

\section{Alternative packages}
\label{sec:alternatives}

The Comprehensive R Archive Network\footnote{%
  \url{https://cran.r-project.org}} %
(CRAN) contains a number of packages with features overlapping
\pkg{expint}. I review the similarities and differences here.

The closest package in functionality is \pkg{gsl} \citep{gsl-package}.
This package is an R wrapper for the special functions and quasi
random number generators of the GNU Scientific Library. As such, it
provides access to basically the same C code as used in \pkg{expint}.
Apart from the changes to the GSL code mentioned in
\autoref{sec:implementation}, the main difference between the two
packages is that installation from source of \pkg{gsl} requires that
the GSL be installed on one's system, whereas \pkg{expint} is a
regular, free standing R package.

\pkg{VGAM} \citep{VGAM} is a large, high quality package that provides
functions to compute the exponential integral $\Ei(x)$ for real
values, as well as $e^{-x} \Ei(x)$ and $E_1(x)$ and their derivatives
(up to the third derivative). Functions \code{expint},
\code{expexpint} and \code{expint.E1} are wrappers to the
Netlib\footnote{%
  \url{https://www.netlib.org}} %
FORTRAN subroutines in file \code{ei.f}. \pkg{VGAM} does not provide
an API to its C routines.

The package \pkg{pracma} \citep{pracma} provides a large number of
functions from numerical analysis, linear algebra, numerical
optimization, differential equations and special functions. Its
versions of \code{expint}, \code{expint\_E1}, \code{expint\_Ei} and
\code{gammainc} are entirely written in R with perhaps less focus on
numerical accuracy than the GSL and Netlib implementations. The
functions are not vectorized, and the incomplete gamma function is
supported for $a \geq -1$ only.

The package \pkg{frmqa} had a function \code{gamma\_inc\_err} that
computed the incomplete gamma function using the incomplete Laplace
integral, but it was only valid for $a = j + 0.5$,
$j = 0, 1, 2, \dots$. (The package was removed from CRAN in 2022.)

Package \pkg{zipfR} \citep{zipfR} introduces a set of functions to
compute various quantities related to the gamma and incomplete gamma
functions, but these are essentially wrappers around the base
R functions \code{gamma} and \code{pgamma} with no new
functionalities.


\section{Examples}
\label{sec:examples}

Let us tabulate the values of $E_n(x)$ for $x = 1.275, 10, 12.3$ and
$n = 1, 2, \dots, 10$ as found in examples~4--6 of
\citet[section~5.3]{Abramowitz:1972}.
<<echo=TRUE>>=
x <- c(1.275, 10, 12.3)
n <- 1:10
structure(t(outer(x, n, expint)),
          dimnames = list(paste("n =", n),
                          paste("x =", x)))
@

We may also tabulate the values of $\Gamma(a, x)$ for
$a = -1.5, -1, -0.5, 1$ and $x = 1, 2, \dots, 10$.
<<echo=TRUE>>=
a <- c(-1.5, -1, -0.5, 1)
x <- 1:10
structure(t(outer(a, x, gammainc)),
          dimnames = list(paste("x =", x),
                          paste("a =", a)))
@


\section{Acknowledgments}

I built on the source code of R and many of the packages cited in
this manual to create \pkg{expint}, so the R Core Team and the
package developers deserve credit. I also extend my thanks to Dirk
Eddelbuettel who was of great help in putting together the package
API, through both his posts in online forums and private
communication. Joshua Ulrich provided a fix to the API infrastructure
to avoid duplicated symbols that was implemented in version 0.1-6 of
the package.

\appendix

\section{Additional details on the computation of the incomplete gamma
  function}
\label{sec:computation_gamma_inc}

\href{https://gitlab.com/vigou3/expint/-/issues/2}{Issue \#2} raised
by Geoffrey Poole highlights that the function \code{gammainc} in
versions of \pkg{expint} prior to 0.2-0 returned inconsistent results
for small negative values of $a$ and a small value of $x$.
\autoref{fig:gammainc_vs_incgam} illustrates the problem by comparing
the behaviour of $\Gamma(a, 10^{-5})$ around $a = -0.5$ between
\code{gammainc} of \pkg{expint} (again, prior to 0.2-0) and
\code{incgam} of \pkg{pracma}.

\begin{figure}[t]
  \centering
<<echo=FALSE, fig=TRUE>>=
a <- seq(-0.501, -0.499, length.out = 1000)
x <- 1e-5
if (requireNamespace("gsl"))
    gGamma <- gsl::gamma_inc(a, x)
if (requireNamespace("pracma"))
    pGamma <- sapply(a, pracma::incgam, x = x)
plot(a, gGamma, type = "l", ylab = expression(Gamma(a, x)),
     main = "Versions of package expint prior to 0.2-0")
lines(a, pGamma, lty = 2, lwd = 2, col = "orange")
legend(-0.5002, 635,
       c("expint::gammainc(a, 1e-05)", "pracma::incgam(1e-05, a)"),
       lty = c(1, 2), lwd = c(1, 2), col = c("black", "orange"))
@
  \caption{Incomplete gamma function for small negative values of $a$
    and a small value of $x$}
  \label{fig:gammainc_vs_incgam}
\end{figure}

The function \code{gamma\_inc} for the package \pkg{gsl} shows the
same defect\footnote{%
  The solid line in Figure \ref*{fig:gammainc_vs_incgam} is actually
  traced using \code{gamma\_inc}.}, %
indicating that the problem must lie in the GSL code. Indeed, the GSL
routine for the incomplete gamma function treats specially the case
$-0.5 < a < 0$, and reverts to the recursion
\eqref{eq:gammainc_recursion} only for $a \leq -0.5$.

For $-0.5 < a < 0$, the GSL computes the incomplete gamma function
using the relation
\begin{equation*}
  \Gamma(a, x) = \Gamma(a) Q(a, x),
\end{equation*}
with $Q(a, x)$ defined as follows:
\begin{align*}
  Q(a, x)
  &= 1 - P(a, x) \\
  &= 1 - \frac{\gamma(a, x)}{\Gamma(a)},
\end{align*}
with
\begin{equation*}
  \gamma(a, x) = P(a, x) \Gamma(a) = \int_0^x t^{a - 1} e^{-t}\, dt.
\end{equation*}
The GSL routine \code{gamma\_inc\_Q\_series} carries the computation
of $Q(a, x)$ using a series expansion. The code is not obvious and
requires a fair share of reverse engineering. Hence, I document my
findings here.

First, a series expansion for $\gamma(a, x)$ is \citep[section
6.5.33]{Abramowitz:1972}:
\begin{equation}
  \gamma(a, x) =
  \sum_{n = 0}^\infty \frac{(-1)^n}{a + n} \frac{x^{a + n}}{n!}.
\end{equation}
We can rewrite this expansion as
\begin{align}
  \gamma(a, x)
  &= \frac{x^a}{a} -
    \sum_{n = 1}^\infty \frac{(-1)^{n - 1}}{a + n} \frac{x^{a + n}}{n!}
    \notag \\
  &= \frac{x^a}{a} -
    \frac{x^{a + 1}}{a + 1}
    \sum_{n = 1}^\infty (-1)^{n - 1}
    \left( \frac{a + 1}{a + n} \right) \frac{x^{n - 1}}{n!}.
    \label{eq:expansion_small_gamma}
\end{align}

Second, using the fact that $\Gamma(a + 1) = a \Gamma(a)$, we can
rewrite $Q(a, x)$ as follows:
\begin{align*}
  Q(a, x)
  &= 1 - \frac{x^a}{\Gamma(a + 1)} + \frac{x^a}{\Gamma(a + 1)}
    - \frac{\gamma(a, x)}{\Gamma(a)} \\
  &= 1 - \frac{x^a}{\Gamma(a + 1)} +
    \frac{1}{\Gamma(a)}
    \left(
      \frac{x^a}{a} - \gamma(a, x)
    \right) \\
  &= 1 - \frac{x^a}{\Gamma(a + 1)} +
    \frac{x^a}{\Gamma(a + 1)}
    \frac{\Gamma(a + 1)}{\Gamma(a)}
    \frac{x}{x^{a + 1}}
    \left(
      \frac{x^a}{a} - \gamma(a, x)
    \right) \\
  &= 1 - \frac{x^a}{\Gamma(a + 1)} +
    \frac{x^a}{\Gamma(a + 1)}
    \left(
      \frac{a}{a + 1}
    \right)
    x
    \left(
      \frac{a + 1}{x^{a + 1}}
    \right)
    \left(
      \frac{x^a}{a} - \gamma(a, x)
    \right).
\end{align*}
If we define $\phi(a) = 1 - x^a/\Gamma(a + 1)$ and we simplify the last
term above using \eqref{eq:expansion_small_gamma}, we finally obtain:
\begin{equation}
  \label{eq:gsl_q_series}
  Q(a, x) = \phi(a) + (1 - \phi(a))
  \left(
    \frac{a}{a + 1}
  \right)
  x
  \sum_{n = 1}^\infty (-1)^{n - 1}
  \left( \frac{a + 1}{a + n} \right) \frac{x^{n - 1}}{n!}.
\end{equation}
This is the expression used to compute $Q(a, x)$.

The only remaining element is the computation of $\phi(a)$. For this,
the routine \code{gamma\_inc\_Q\_series} uses the product of the
Taylor series expansion of $x^a$ around $a = 0$,
\begin{equation*}
  x^a = \sum_{n = 0}^\infty \frac{(a \ln x)^n}{n!},
\end{equation*}
and the following Taylor series expansion of the reciprocal gamma
function \citep{Mathworld:Gamma}:
\begin{equation*}
  \frac{1}{\Gamma(a + 1)} = 1 + \gamma a +
  \left( \frac{\gamma^2}{2} - \frac{\pi^2}{12} \right) a^2 +
  \left(
    \frac{\gamma^3}{6} - \frac{\gamma\pi^2}{12} + \frac{\zeta(3)}{3}
  \right) a^3 + \cdots,
\end{equation*}
where $\gamma$ is the Euler constant, and $\zeta$ is the Riemann zeta
function \citep[chapters 6 and 23]{Abramowitz:1972}. Putting all these
pieces together, we obtain:
\begin{align*}
  \phi(a)
  &= 1 - \frac{x^a}{\Gamma(a + 1)} \\
  &= - (\ln x + \gamma) a +
    \left\{
      \frac{\pi}{12} - (\ln x + \gamma)^2
    \right\} a^2 \\
  &\phantom{=} +
    \left\{
      (\ln x + \gamma)
      \left(
        \frac{\pi}{12} - \frac{(\ln x + \gamma)^2}{6}
      \right) -
      \frac{\zeta(3)}{3}
    \right\} a^3 + \cdots.
\end{align*}
One will recognize above the coefficients \code{c1}, \code{c2} and
\code{c3} of \code{term1} in the routine \code{gamma\_inc\_Q\_series}.
It is not clear where the numerical values in the other seven
coefficients come from. Here endeth reverse engineering.

The fix introduced in \pkg{expint} 0.2-0 follows the same strategy as
\pkg{pracma}: just use the recursion \eqref{eq:gammainc_recursion}
also for $-0.5 < a < 0$, and rely on the accuracy of \code{pgamma} for
small values of $a$ to yield the correct result.
\autoref{fig:gammainc_vs_incgam_fixed} shows that the results from
\pkg{expint} and \pkg{pracma} now match.

One last implementation detail: the original GSL code uses a loop to
compute \eqref{eq:gammainc_recursion} as many times as necessary for
$a < -0.5$. Simply extending usage of the loop to $a < 0$ does not
work due to rounding errors in computations involving values of $a$ in
$(-0.5, 0)$. Therefore, \pkg{expint} keeps as a special case the
single application of the recursive relation for $a$ in the latter
range.

\begin{figure}[t]
  \centering
<<echo=FALSE, fig=TRUE>>=
a <- seq(-0.501, -0.499, length.out = 1000)
x <- 1e-5
eGamma <- gammainc(a, x)
if (requireNamespace("pracma"))
    pGamma <- sapply(a, pracma::incgam, x = x)
plot(a, eGamma, type = "l", ylab = expression(Gamma(a, x)),
     main = "Versions of package expint from 0.2-0")
lines(a, pGamma, lty = 2, lwd = 2, col = "orange")
legend(-0.5002, 635,
       c("expint::gammainc(a, 1e-05)", "pracma::incgam(1e-05, a)"),
       lty = c(1, 2), lwd = c(1, 2), col = c("black", "orange"))
@
  \caption{Fixed incomplete gamma function}
  \label{fig:gammainc_vs_incgam_fixed}
\end{figure}

\bibliography{expint}

\end{document}