\documentclass[x11names,english]{article}
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
  \usepackage{amsmath}
  \usepackage[round]{natbib}
  \usepackage{babel}
  \usepackage{microtype}
  \usepackage[scaled=0.92]{helvet}
  \usepackage[sc]{mathpazo}
  \usepackage[noae,inconsolata]{Sweave}
  \usepackage{framed}

  %\VignetteIndexEntry{expint user manual}
  %\VignettePackage{expint}
  %\SweaveUTF8

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

  %% Colors
  \usepackage{xcolor}
  \definecolor{link}{rgb}{0,0.4,0.6}             % internal links
  \definecolor{url}{rgb}{0.6,0,0}                % external links
  \definecolor{citation}{rgb}{0,0.5,0}           % citations
  \definecolor{codebg}{named}{LightYellow1}      % R code background

  %% Hyperlinks
  \usepackage{hyperref}
  \hypersetup{%
    pdfauthor={Vincent Goulet},
    colorlinks = {true},
    linktocpage = {true},
    urlcolor = {url},
    linkcolor = {link},
    citecolor = {citation},
    pdfpagemode = {UseOutlines},
    pdfstartview = {Fit},
    bookmarksopen = {true},
    bookmarksnumbered = {true},
    bookmarksdepth = {subsubsection}}

  %% Sweave Sinput and Soutput environments reinitialized to remove
  %% default configuration. Space between input and output blocks also
  %% reduced.
  \DefineVerbatimEnvironment{Sinput}{Verbatim}{}
  \DefineVerbatimEnvironment{Soutput}{Verbatim}{}
  \fvset{listparameters={\setlength{\topsep}{0pt}}}

  %% Environment Schunk redefined as an hybrid of environments
  %% snugshade* and leftbar of framed.sty.
  \makeatletter
  \renewenvironment{Schunk}{%
    \setlength{\topsep}{1pt}
    \def\FrameCommand##1{\hskip\@totalleftmargin
       \vrule width 2pt\colorbox{codebg}{\hspace{3pt}##1}%
      % There is no \@totalrightmargin, so:
      \hskip-\linewidth \hskip-\@totalleftmargin \hskip\columnwidth}%
    \MakeFramed {\advance\hsize-\width
      \@totalleftmargin\z@ \linewidth\hsize
      \advance\labelsep\fboxsep
      \@setminipage}%
  }{\par\unskip\@minipagefalse\endMakeFramed}
  \makeatother

  %% Additional commands similar to those of RJournal.sty.
  \newcommand{\pkg}[1]{\textbf{#1}}
  \newcommand{\code}[1]{\texttt{#1}}
  \newcommand{\samp}[1]{{`\normalfont\texttt{#1}'}}
  \newcommand{\file}[1]{{`\normalfont\textsf{#1}'}}

  \newcommand{\Ei}{\operatorname{Ei}}

  \bibliographystyle{plainnat}

<<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 that intends to fill a gap in R's
support for mathematical functions by providing 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}.

Package \pkg{expint} started its life in version~2.0-0 of package
\pkg{actuar} \citep{actuar} where we 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
relation
\begin{equation*}
  \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}

Package \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. 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)
@

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

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

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, 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}

Package \pkg{expint} makes these routines available to other packages
through declarations in the header file \file{include/expintAPI.h} in
the package installation directory. The developer of some package
\pkg{pkg} who wants to use a routine --- say \code{expint\_E1} --- in
her code should proceed as follows.
\begin{enumerate}
\item Add 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 that will call
  \code{expint\_E1}, say \code{pkg\_expint\_E1} to avoid any name
  clash, in \file{src/init.c} as follows:
\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 \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 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 package \pkg{zoo} \citep{zoo}. 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 no 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 we conducted
proved the approach we 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, we made the following changes
to the original GSL code:
\begin{enumerate}
\item routines now compute a single value and return their result by
  value;
\item accordingly, calculation of the approximation error was dropped
  in all routines;
\item most importantly, \code{gamma\_inc} does not compute
  $\Gamma(a, x)$ for $a > 0$ with \eqref{eq:gammainc_vs_pgamma} using
  the GSL routines, but rather using routines \code{gammafn} and
  \code{pgamma} part of the R API.
\end{enumerate}
The following illustrates the last point.
<<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
@

\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}. We 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 \pkg{gsl} requires that the GSL be installed on one's
system, whereas \pkg{expint} is a regular, free standing R
package.

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{http://www.netlib.org}} %
FORTRAN subroutines in file \code{ei.f}. \pkg{VGAM} does
not provide an API to its C routines.

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 version of \code{gammainc} only supports positive
values of $a$.

Package \pkg{frmqa} \citep{frmqa} has a function
\code{gamma\_inc\_err} that computes the incomplete gamma function
using the incomplete Laplace integral, but it is only valid for
$a = j + \frac{1}{2}$, $j = 0, 1, 2, \dots$.

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}

We 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(n, paste("x =", x)))
@

We 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(x, paste("a =", a)))
@


\section{Acknowledgments}

We built on the source code of R and many of the packages cited in
this paper to create \pkg{expint}, so the R Core Team and the package
developers deserve credit. We also extend our 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.


\bibliography{expint}

\end{document}