\section{Numerical integration of polynomial spirals}
\label{sec-numint}

Almost all calculations involving polynomial spiral curves can be
reduced to an evaluation of a single definite integral,

\begin{equation}
\label{spiro-int}
\mbox{\textit{spiro}}(k_0, k_1, k_2, k_3) = \int_{-0.5}^{0.5}
e^{i(k_0 s + \frac{1}{2}k_1s^2 + \frac{1}{6}k_2s^3 +
      \frac{1}{24}k_3s^4)} ds\:.
\end{equation}

Here, $k_0 ... k_3$ are the curvature and its first three derivatives
at the center of a curve segment of arc length 1. More generally (for
an arbitrary degree polynomial),

\begin{equation}
\label{spiro-int2}
\mbox{\textit{spiro}}(\textbf{k}) = \int_{-0.5}^{0.5}
e^{i\sum_{0 \leq j}\frac{1}{(j + 1)!}k_js^{j + 1}} ds\:.
\end{equation}

Fixing the limits to $\pm 0.5$ simplifies the math and (as will be
seen) creates symmetries which can be exploited to improve numerical
accuracy, but does not lose generality. An integral over any limits
can be expressed in terms of the spiro integral, as follows:

\begin{equation}
\label{spiro-arbitrary-limits}
\begin{split}
\int_{s_0}^{s_1}
e^{i(k_0 s + \frac{1}{2}k_1s^2 + \frac{1}{6}k_2s^3 +
      \frac{1}{24}k_3s^4)} ds = \\[2mm]
(s_1 - s_0)e^{i(k_0s + \frac{1}{2}k_1s^2 + \frac{1}{6}k_2s^3 +
  \frac{1}{24}k_3s^4)} \\
\mbox{\textit{spiro}}( & (s_1 - s_0)(k_0 + k_1s + \tfrac{1}{2}k_2s^2 + \tfrac{1}{6}k_3s^3), \\
& (s_1 - s_0)^2(k_1 + k_2s + \tfrac{1}{2}k_3s^2), \\
& (s_1 - s_0)^3(k_2 + k_3s), \\
& (s_1 - s_0)^4k_3),\  \\[2mm]
\mbox{where}\ s = \frac{s_0 + s_1}{2}\:.
\end{split}
\end{equation}

Note that, as a special case of Equation \ref{spiro-arbitrary-limits},
the classic Fresnel integrals also reduce readily to this integral.

\begin{equation}
\int_0^s e^{it^2} dt = se^{i\frac{s^2}{8}}
\textit{spiro}(\tfrac{s}{2}, s^2, 0, 0)\:.
\end{equation}

The spiro integral can be thought of as analogous to the sine and
cosine functions of the generalized polynomial spiral. This connection
is clear when considering a circular arc.

\begin{equation}
\frac{\sin x}{x} = \textit{spiro}(2s, 0, 0, 0)\:.
\end{equation}

The absolute value of the spiro integral is the ratio of chord length to
arc length of the curve segment (note that expressing this quantity as
a ratio is invariant to scaling, and it is obviously invariant to
rotation and translation as well). For a given curve segment with
chord length $\Delta x$, the curvature and
its derivatives along the curve (with arc length $s$ normalized to the
range $[-0.5, 0.5]$) are given by a simple formula,

\begin{equation}
\kappa^{(j)}(s) = \Big(\frac{|\textit{spiro}(\textbf{k})|}{\Delta x}\Big)^{j+1}\sum_{0 \leq
  i} \frac{k_{i + j}}{(i + 1)!}s^i\:.
\end{equation}

\subsection{Polynomial approximation of the integral}

There are a number of approaches to numerically computing an integral
such as Equation \ref{spiro-int}. However, a straightforward numerical
integration technique such as Runge--Kutte would converge only
relatively slowly, requiring significant computation time to achieve
high accuracy. The approach of this section is to compute polynomial
approximations of the integral, keeping all terms up to a certain
order and discarding the rest.

A key insight is that the angle $\theta(s)$ is a straightforward
polynomial function of $s$. Thus, to integrate $e^{i\theta(s)}$, take
the Taylor's series expansion of $e^x$, substitute the polynomial
$i\theta(s)$ for $x$, and expand symbolically, keeping terms only up
to the given order.

\begin{equation}
\label{theta-closedform}
\theta(s) = k_0 s + \frac{k_1}{2}s^2 + \frac{k_2}{6}s^3 +
      \frac{k_3}{24}s^4\:.
\end{equation}

Below is a worked example retaining terms for $e^{ix}$ up to
quadratic, then integrating. The approximation is

\begin{equation}
e^{ix} = 1 + ix - \frac{1}{2}x^2 \hspace{10mm} + O(x^3) \:.
\end{equation}

Substituting $x = \theta(s)$, we get

\begin{equation}
e^{i\theta(s)} = 1 + ik_0 s - \frac{i k_1 - k_0^2}{2}
  s^2 \hspace{10mm} + O(s^3)\:.
\end{equation}

Integrating, we get

\begin{equation}
\int e^{i\theta(s)}ds = f(s) = s + \frac{ik_0}{2}s^2 - \frac{i k_1 - k_0^2}{6}
  s^3 \hspace{10mm} + O(s^4)\:.
\end{equation}

Evaluating $f(0.5) - f(-0.5)$, we find that all even
polynomial terms (including $s^4$) drop out:

\begin{equation}
f(0.5) - f(-0.5) \approx 1 + \frac{ik_1 - k_0^2}{24}\:.
\end{equation}

Thus, as promised, fixing the limits of integration to $\pm 0.5$
provides a symmetry that both decreases the number of terms to compute
and also increases the accuracy by one degree.

Retaining all terms up to quartic yields an approximation not much
more complicated. Note that this approximation actually uses all four
parameters.

\begin{equation}
\Delta x \approx 1 + \frac{ik_1 - k_0^2}{24} + \frac{ik_3 - 6ik_0^2k_1
    + k_0^4-4k_0k_2 - 3k_1^2}{1920}\:.
\end{equation}

\subsection{Higher-order polynomial approximations}

Higher-order polynomials offer a good way to increase the accuracy of
the approximation at modest cost. However, manually working out the
coefficients becomes quite tedious. This section describes a technique
for automatically generating code for an arbitrary polynomial degree.

For functions of a single variable, the Taylor's series expansion
gives the coefficients of the approximating polynomial simply and
unambiguously (it is, of course, possible to use refinements such as
Chebysheff polynomials to redistribute the errors more evenly, but the
underlying concept remains). Here, the task is to approximate a
function of four parameters. It is possible to generalize the Taylor's
series expansion to functions of more than one parameter, but a
simplistic approach would result in a blowup of coefficients; for a
function of four parameters, the number of coefficients grows
quadratically with the order of the polynomial.

A more sophisticated approach is represent each successive $x^j$ as a
vector of polynomial coefficients. Because of the symmetries of the
equation (many odd terms will drop out), for $j \ge 3$, it is most
efficient to compute $x^j$ by combining the vectors for $x^{j-2}$ and
$x^2$. Then, these polynomials are substituted into the Taylor's
series expansion for $e^x$. A fully worked out example, for order 8,
is given below.

\begin{equation}
\label{order-8-spiro}
\begin{array}{l}
t_{11} = k_0 \\[-1mm]
t_{12} = \frac{k_1}{2} \\[-1mm]
t_{13} = \frac{k_2}{6} \\[-1mm]
t_{14} = \frac{k_3}{24} \\[-1mm]
t_{22} = t_{11}^2 \\[-1mm]
t_{23} = 2 t_{11} t_{12} \\[-1mm]
t_{24} = 2 t_{11} t_{13} + t_{12}^2 \\[-1mm]
t_{25} = 2 (t_{11} t_{14} + t_{12} t_{13}) \\[-1mm]
t_{26} = 2 t_{12} t_{14} + t_{13}^2 \\[-1mm]
t_{27} = 2 t_{13} t_{14} \\[-1mm]
t_{28} = t_{14}^2 \\[-1mm]
t_{34} = t_{22} t_{12} + t_{23} t_{11} \\[-1mm]
t_{36} = t_{22} t_{14} + t_{23} t_{13} + t_{24} t_{12} + t_{25} t_{11} \\[-1mm]
t_{38} = t_{24} t_{14} + t_{25} t_{13} + t_{26} t_{12} + t_{27} t_{11} \\[-1mm]
t_{44} = t_{22}^2 \\[-1mm]
t_{45} = 2 t_{22} t_{23} \\[-1mm]
t_{46} = 2 t_{22} t_{24} + t_{23}^2 \\[-1mm]
t_{47} = 2 (t_{22} t_{25} + t_{23} t_{24}) \\[-1mm]
t_{48} = 2 (t_{22} t_{26} + t_{23} t_{25}) + t_{24}^2 \\[-1mm]
t_{56} = t_{44} t_{12} + t_{45} t_{11} \\[-1mm]
t_{58} = t_{44} t_{14} + t_{45} t_{13} + t_{46} t_{12} + t_{47} t_{11} \\[-1mm]
t_{66} = t_{44} t_{22} \\[-1mm]
t_{67} = t_{44} t_{23} + t_{45} t_{22} \\[-1mm]
t_{68} = t_{44} t_{24} + t_{45} t_{23} + t_{46} t_{22} \\[-1mm]
t_{78} = t_{66} t_{12} + t_{67} t_{11} \\[-1mm]
t_{88} = t_{66} t_{22} \\[-1mm]
u = 1 - \frac{t_{22}}{24} - \frac{t_{24}}{160} - \frac{t_{26}}{896} - \frac{t_{28}}{4608} + \frac{t_{44}}{1920} + \frac{t_{46}}{10752} + \frac{t_{48}}{55296} - \frac{t_{66}}{322560} - \frac{t_{68}}{1658880} + \frac{t_{88}}{92897280} \\[-1mm]
v = \frac{t_{12}}{12} + \frac{t_{14}}{80} - \frac{t_{34}}{480} - \frac{t_{36}}{2688} - \frac{t_{38}}{13824} + \frac{t_{56}}{53760} + \frac{t_{58}}{276480} - \frac{t_{78}}{11612160} \\[-1mm]
\end{array}
\end{equation}

Obviously, for higher degrees, the formulas are too complex to be
worked out by hand. The accompanying software distribution includes a
utility, \texttt{numintsynth.py}, that synthesizes C code for an
arbitrary degree. Asymptotically, the number of lines of code grows as
$O(n^2)$, so for any reasonably low degree, the code is quite
concise. Additionally, the structure of the dependency graph permits a
high level of instruction-level parallelism, making the computation
very efficient on typical CPUs.

\subsection{Hybrid approach}

% spiro-int is not in makefile. to make spiro-int.ps:
% sh golf/curves/run-spiro-int
\begin{figure*}
\label{fig-spiro-timings}
\begin{center}
\includegraphics[scale=.6]{figs/spiro-int}
\caption{Precision vs. timing for different integration strategies.}
\end{center}
\end{figure*}

The higher order the polynomial, the more accurate the
result. However, code complexity grows significantly, and for small
angles, accurate results are possible with a low-order polynomial.

A hybrid approach is to implement a polynomial of fixed order, and use
it for parameter values where the error is within tolerance. For
parameter values outside this tolerance, subdivide. The range $[-0.5,
  0.5]$ is split into $n$ equal subdivisions, each is evaluated using
Equation \ref{spiro-arbitrary-limits}, and the result is the sum of
all such intervals.

The extreme case of an order-2 polynomial is equivalent to Euler
integration. The code is simple, but convergence is expected to be
poor.

The optimum order of polynomial is determined empirically. All even
orders of polynomial between $2$ and $16$ were measured, and the
accuracy plotted against CPU time. The results are shown in Figure
\ref{fig-spiro-timings}, measured on a MacBook Pro with a 2.33 GHz Core
2 Duo processor. Each successive datapoint doubles the total number of
subdivisions. This particular plot is for evaluating $\textit{spiro}(1,
2, 3, 4)$, but other combinations of arguments all give similar
results.

The $n=1$ case is relatively faster than $n>1$, because no sine and
cosine functions are needed for the subdivision. As can be seen, low
order polynomials require a large number of subdivisions to provide
accurate results. However, there is diminishing return as the order of
the polynomial grows -- while the number of intervals needed is
decreased, the time needed to compute each one increases. As can be
seen from the graph, orders greater than 12 don't improve overall
performance significantly.

An empirically determined estimate for the error of the order 12 spiro
approximation is:

\begin{equation}
|\textit{spiro}_{12}(k_0, k_1, k_2, k_3) - \textit{spiro}(k_0,
k_1, k_2, k_3)|
\approx (0.006|k_0|^2 + 0.03|k_1| + 0.03|k_2|^{\frac{2}{3}} +
0.025|k_3|^{\frac{1}{2}})^6
\end{equation}

With $n$ subdivisions, the error scales approximately as
$n^{-12}$. Thus, it should be straightforward to compute a value for
any given precision.