  We provide in this work an algorithm for approximating a very broad class of
symmetric Toeplitz matrices to machine precision in $\mathcal{O}(n \log n)$
time with applications to fitting time series models. In particular, for a
symmetric Toeplitz matrix $\mathbf{\Sigma}$ with values $\mathbf{\Sigma}_{j,k}
= h_{|j-k|} = \int_{-1/2}^{1/2} e^{2 \pi i |j-k| \omega} S(\omega) \mathrm{d}
\omega$ where $S(\omega)$ is piecewise smooth, we give an approximation
$\mathbf{\mathcal{F}} \mathbf{\Sigma} \mathbf{\mathcal{F}}^H \approx \mathbf{D}
+ \mathbf{U} \mathbf{V}^H$, where $\mathbf{\mathcal{F}}$ is the DFT matrix,
$\mathbf{D}$ is diagonal, and the matrices $\mathbf{U}$ and $\mathbf{V}$ are in
$\mathbb{C}^{n \times r}$ with $r \ll n$. Studying these matrices in the
context of time series, we offer a theoretical explanation of this structure
and connect it to existing spectral-domain approximation frameworks. We then
give a complete discussion of the numerical method for assembling the
approximation and demonstrate its efficiency for improving Whittle-type
likelihood approximations, including dramatic examples where a correction of
rank $r = 2$ to the standard Whittle approximation increases the accuracy from
$3$ to $14$ digits for a matrix $\mathbf{\Sigma} \in \mathbb{R}^{10^5 \times
10^5}$. The method and analysis of this work applies well beyond time series
analysis, providing an algorithm for extremely accurate direct solves with a
wide variety of symmetric Toeplitz matrices. The analysis employed here largely
depends on asymptotic expansions of oscillatory integrals, and also provides a
new perspective on when existing spectral-domain approximation methods for
Gaussian log-likelihoods can be particularly problematic.
