Transmutagen implements the Remez algorithm to compute the CRAM approximation.

All computations are performed using SymPy~\cite{10.7717/peerj-cs.103}

symbolic expressions with arbitrary precision floating point numbers, which

are backed by the mpmath library~\cite{ationneeded}.

To compute the CRAM approximation, the interval $[0, \infty)$ is first

translated to $[-1, 1)$ via the transformation $t\mapsto c\frac{t+1}{t-1}$.

\todo{cite Carpenter paper here.} \todo{discuss this {\it after} discussing

finding the roots.} This is done so that the maxima of the error at each

stage are all on a bounded interval (see below). The constant $c$ is chosen so

that these are distributed sufficiently evenly across the interval $[-1, 1)$.

If they are clustered too much, the numeric root finding algorithm may fail to

find them all. We found that the value $c=k/\phi$, where $\phi=1.618\ldots$ is

the golden ratio, works well. The choice of $c$ does not affect the final

result on $[0, \infty)$, but deviations from this heuristic resulted in issues

in the root finding stage for large degrees.

To help visualize this, Figure~\ref{fig:cram-plot} shows the error in the

shifted approximation on $[-1, 1)$ and the final translated approximation on

$[0, \infty)$, for degree 14.

\begin{figure}[!ht]

\resizebox{\textwidth}{!}{\input{cram-plot.pgf}}

\caption{Left: plot of $\hat{r}_{14, 14}\left(c\frac{t+1}{t-1}\right) -

e^{-c\frac{t+1}{t-1}}$; right: plot of $\hat{r}_{14, 14}(t) -

e^{-t}$ ($c=\frac{14}{\phi}\approx 8.652$).}

\label{fig:cram-plot}

\end{figure}

In transmutagen's implementation of the Remez algorithm, the expressions

\begin{equation}

\end{equation}

\begin{equation}

\end{equation}

and

\begin{equation}

\end{equation}

are represented symbolically as SymPy expressions. The expression $E$ is then

evaluated in $t$ at $2(k + 1)$ points in the interval $(-1, 1)$ for

$i=0\ldots 2k+1$. For the first iteration of the algorithm, any set of initial

points in $(-1, 1)$ can be used. By default, transmutagen uses the Chebyshev

nodes~\cite{ationneeded}, that is, the roots of $T_{2(k +1)}(x)$, which always

lie in the interval $(-1, 1)$.\footnote{$T_n(x)$ is the $n$th Chebyshev

polynomial of the first kind, defined as $T_n(\cos(x)) = \cos(nx)$, e.g.,

$T_2(x) = 2x^2 - 1$. The roots of $T_n(x)$ are

$\cos{\left (\frac{2k - 1}{n}\frac{\pi}{2} \right )},\,k=1\ldots n$.} It can

also optionally use a random set of points. We found that convergence can take

as much as 50\% more steps when random initial points are used instead of the

Chebyshev nodes. This evaluation of $E$ results in a

nonlinear system of $2(k+1)$ equations in $2(k+1)$ variables

($p_0,\ldots,p_k,q_1,\ldots,q_k,\epsilon$). These are solved using the

\texttt{sympy.nsolve} function, which internally uses \texttt{mpmath.findroot}

to apply Newton's method~\cite{ationneeded} with a symbolically computed

Jacobian.

\label{sec:nsolve-bisection}

The solution to this system is then substituted into $D$, and the critical

points of this function on the interval $[-1, 1)$ are then found. This is done

by taking the symbolic derivative of $D$, splitting the interval into

subintervals, and finding the roots on each subinterval using the bisection

algorithm, via \texttt{nsolve}. Let $\{z_i\}$ be the set of critical points on

$[-1, 1)$, including the end-points $D|_{t=-1}$ and

$\lim_{t\to 1} D=-r|_{t=1}$. There are $2(k+1)$ points $\{z_i\}$. These points

$\{z_i\}$ are used as the set of initial points for the next iteration, and

the algorithm iterates as such until convergence is reached.

By the equioscillation theorem~\cite{ationneeded}, the approximation

$\hat{r}_{k, k}$ is minimal for a given order $k$ when the critical points of

$\hat{r}_{k, k}(t) - e^{-t}$ oscillate in sign and have equal absolute value.

For the iterates of the algorithm, the points $z_i$ are these critical points

of the current approximation, so convergence is detected by looking at

$\varepsilon_N = \max{|z_i|} - \min{|z_i|}$, for the $N$th iteration of the

algorithm. Convergence occurs roughly when $\varepsilon_N$ is near $10^{-d}$,

where $d$ decimal digits are used in the calculations. However, the true

minimal value of $\varepsilon_N$ depends on both $k$ and $d$. We found a

robust heuristic to be to iterate until $\log_{10}{(\varepsilon_N)}$ is near

$-d$, then stop iterating when the values of $\varepsilon_N$ become

log-convex. See Figure~\ref{fig:convergence-14-1000} for an example of the

convergence of $\varepsilon_N$ for degree 14, 1000 digits of precision.

\begin{figure}[!ht]

\resizebox{0.9\textwidth}{!}{\input{convergence-14-1000.pgf}}

\caption{Convergence of the Remez algorithm for degree 14, 1000 digits of

precision.}

\label{fig:convergence-14-1000}

\end{figure}

Once the algorithm converges, the rational function is translated back to the

interval $[0, \infty)$ via the inverse transformation

$t\mapsto \frac{t - c}{t + c}$. This must be done with care to avoid losing

precision. For a polynomial $f$ and rational function $p/q$, the SymPy method

\texttt{Poly(f).\allowbreak{}transform(p, q)} efficiently computes

$q^nf\left(p/q\right)$ without losing precision. The resulting rational

function is normalized so that the constant term in the denominator is 1,

which matches the form used in the literature.~\cite{ationneeded}

The Remez algorithm is outlined in Algorithm~\ref{alg:remez-pseudocode} as pseudocode.

\begin{algorithm}

\caption{The Remez algorithm for the CRAM approximation of $e^{-t}$ on

$[0, \infty)$ of degree $k,k$.}\label{alg:remez-pseudocode}

\begin{algorithmic}[1]

\REQUIRE $k$, the order of the approximation

\STATE $t, p_0, \ldots, p_k, q_1, \ldots, q_k, \epsilon$ are symbolic variables

\STATE $r \leftarrow \frac{p_kt^k + \cdots + p_1t + p_0}{q_kt^k + \cdots +

\STATE $c \leftarrow \frac{k}{\phi}$

\STATE $D \leftarrow e^{c\frac{t+1}{t-1}} - r$

\STATE $\{z_i\} \leftarrow$ Chebyshev nodes of order $2(k+1)$

\STATE $N \leftarrow 1$

\REPEAT

\FOR {$i=1$ \TO $2(k+1)$}

\STATE $E_i \leftarrow D|_{t=z_i} + (-1)^i\epsilon$

\ENDFOR

\STATE $sol \leftarrow$ Solve $\{E_i\}$ for $(p_0,\ldots,p_k,q_1,\ldots,q_k,\epsilon)$

\STATE $z_0 \leftarrow D|_{sol,t=-1}$

\STATE $z_i \leftarrow$ $i^\mathrm{th}$ critical point of $D|_{sol}$ on

$(-1, 1)$ \COMMENT{$i=1\ldots 2n$}

\STATE $z_{2(k + 1)} \leftarrow -r|_{sol,t=1}$ \COMMENT{$\lim_{t\to 1} D|_{sol}$}

\STATE $\varepsilon_N \leftarrow \max{|z_i|} - \min{|z_i|}$

\STATE $N \leftarrow N + 1$

\UNTIL {convergence condition: $\log_{10}(\varepsilon_N) \approx -d$ \AND

$\varepsilon_N$ is log-convex}

\STATE $r_\mathrm{final}=r|_{t=\frac{t - c}{t + c}}$ \COMMENT{translate

$[-1, 1)$ back to $[0, \infty)$ and normalize $q_0=1$}

\ENSURE $r_\mathrm{final}$

\end{algorithmic}

\end{algorithm}