pusa-comparison.tex

In addition, we have compared our partial fraction digits to those published
in~\cite{pusa2012correction} for degrees 14 and 16. However, many
discrepancies were found. These discrepancies are summarized in
Tables~\ref{table:pusa-degree-14} and~\ref{table:pusa-degree-16}.

Two values, $\mathrm{re}(\alpha_5)$ and $\mathrm{re}(\alpha_7)$ of degree 16,
differ within \texttt{double} machine floating point precision (each within
one-bit). The values cause the approximations to differ by as much as
$6\times10^{-19}$ (see Figure~\ref{fig:pusa-differences}).

\setul{}{1pt}

\begin{table}[h!]
\centering
\begin{tabular}{ r b{1.8in} b{1.8in} }
Coefficient & Values from~\cite{pusa2012correction} & Values computed by Transmutagen \\
\midrule
\input{pusa-table-14.tex}
\bottomrule
\end{tabular}
\caption{Computed partial fraction coefficients for degree 14. Differing
  digits are underlined. \texttt{j} is Python syntax for imaginary numbers
  (e.g., \texttt{2j}$=2i$). Note that~\oldcite{pusa2012correction} only has 19
  digits for $\mathrm{re}(\theta_5)$.}
\label{table:pusa-degree-14}
\end{table}

\begin{table}[h!]
\centering
\begin{tabular}{ r b{1.8in} b{1.8in} }
Coefficient & Values from~\cite{pusa2012correction} & Values computed by Transmutagen \\
\midrule
\input{pusa-table-16.tex}
\bottomrule
\end{tabular}
\caption{Computed partial fraction coefficients for degree 16. Differing
  digits are underlined. \texttt{j} is Python syntax for imaginary numbers
  (e.g., \texttt{2j}$=2i$). Note that $\mathrm{re}(\alpha_5)$ and
  $\mathrm{re}(\alpha_7)$ differ within machine floating point
  precision.}
\label{table:pusa-degree-16}
\end{table}

\begin{figure}[!ht]
\centering
\resizebox{0.9\textwidth}{!}{\input{pusa-differences.pgf}}
\caption{Difference between $\hat{r}_{k,k}$ for the transmutagen computed
  values and the values from~\oldcite{pusa2012correction} for $k=14,16$.}
\label{fig:pusa-differences}
\end{figure}

Let $D(t) = e^{-t} - \hat{r}_{k,k}(t)$. From Equation~\ref{eq:part-frac}, we can see that
$\lim_{t\to\infty}{\left|D(t) \right|} = \alpha_0$. By the
equioscillation theorem, the error of the approximation equioscillates $2k$
times in $(0, \infty)$, and additionally takes on the same value at
$\infty$.\footnote{Recall from Section~\ref{sec:remez-algorithm} that on the
  translated interval $[-1, 1)$, there are $2(k+1)$ critical points of $D$,
  $z_i$, which represent the minimax error of the approximation, including
  $z_0 = D|_{t=-1}$ and $z_{2(k+1)}=\lim_{t\to 1}{D}$. When the approximation
  is translated from $[-1, 1)$ to $[0, \infty)$, the point $z_{2(k + 1)}$
  represents the ``error at $\infty$''. c.f.\ Figure~\ref{fig:cram-plot}.}
What this means is that the computed value for $\alpha_0$ should correspond to
the maximum absolute error of the approximation.

This provides us with a method to test the validity of partial fraction
coefficients. If a value of the difference $D(t)$ exceeds $\alpha_0$ in
absolute value, the coefficients are inconsistent. Conversely, if the critical
points of $D(t)$, which can be found numerically via gradient descent or by
numerically solving the symbolic derivative, are near $\alpha_0$, we can have
a degree of confidence that the coefficients are correct.

Using the subinterval method outlined in Section~\ref{sec:nsolve-bisection}
above, we computed 24 and 27 critical points of $D(t)$ on $[0, 100]$ for
degrees 14 and 16, respectively. These points were then used as the initial
points for \texttt{sympy.\allowbreak{}nsolve}, which internally uses
\texttt{mpmath.findroot}'s secant solver, to solve $\frac{dD}{dt}=0$. This was
done for the partial fraction expressions generated by transmutagen and for
the ones from the coefficients from~\cite{pusa2012correction}. When performing
this analysis, it is important to round the transmutagen partial fraction
coefficients to 20 digits, to give a fair comparison to the coefficients
from~\cite{pusa2012correction}. This is because the differences between the
computed critical values and $\alpha_0$ will themselves have an error on the
order of $10^{-k}$ when using coefficients with $k$ digits.

The results of this analysis are shown in
Figures~\ref{fig:pusa-differences-errors-14}
and~\ref{fig:pusa-differences-errors-16}. For the critical points computed,
the transmutagen partial fraction coefficients produce an error closer to
$\alpha_0$ in each case. For each point, the respective $\alpha_0$ (from
transmutagen or from~\cite{pusa2012correction}) was used, although it should
be noted that the computed $\alpha_0$'s are the same for degree 16, and differ
by only 1 in the least significant digit for degree 14 (see
Tables~\ref{table:pusa-degree-14} and~\ref{table:pusa-degree-16}).

\begin{figure}[!ht]
\centering
\resizebox{0.9\textwidth}{!}{\input{pusa-differences-errors-14.pgf}}
\caption{Difference between $D(t) = e^{-t} - \hat{r}_{14,14}(t)$ and $\alpha_0$ at 24 critical points
$t_0$ of $D(t)$ in $[0, 100]$.}
\label{fig:pusa-differences-errors-14}
\end{figure}

\begin{figure}[!ht]
\centering
\resizebox{0.9\textwidth}{!}{\input{pusa-differences-errors-16.pgf}}
\caption{Difference between $D(t) = e^{-t} - \hat{r}_{16,16}(t)$ and $\alpha_0$ at 27 critical points
$t_0$ of $D(t)$ in $[0, 100]$.}
\label{fig:pusa-differences-errors-16}
\end{figure}