number of timesteps. This is the best behaviour achievable; to do better we would have to move to extended precision or use fewer operations. This fundamental limit is known as Brouwer’s law [@Newcomb1899; @Brouwer1937].

Third, if any parts of the integration algorithm are biased, the errors will be correlated. This leads to a faster long-term energy-error growth than if errors are uncorrelated; it grows linearly with time, i.e. $E_{\rm bias}\sim N \sim t$.

For a given integrator, which of these three error terms dominates depends on the nature of the simulation, the timestep, and the total integration time (number of timesteps).

Improvements {#sec:improvements}
============

The algorithms we describe in Sect. \[sec:background\] have been used successfully for many years. In the following, we show how to significantly improve the speed and accuracy of the algorithms by taking special care in the implementation of several details, many of which are related to finite floating-point precision on modern computers.

For the remainder of this paper, we will assume that we work with a CPU that follows the IEEE 754 standard for floating-point arithmetic. Most importantly, we assume that all floating-point operations follow the *rounding to nearest, ties to even* rule [@IEEE754]. What follows 