Analytically, the traditional way forward is to average over the short-period oscillations in the problem to yield approximate solutions. The great insight of Wisdom and Holman was that, at the same level of approximation, one can [*add*]{} high frequency terms. By judicious choice of these additional frequencies, the perturbations among the planets can be transformed into trivially integrated delta functions. The result is an exceedingly efficient integrator that has proven an indispensable tool for modern studies in celestial mechanics.

In this paper, we present results from a complete reimplementation of the Wisdom-Holman integrator. We show how to speed up the algorithm in several ways and dramatically increase its accuracy. Many of the improvements are related to finite double floating-point precision on modern computers [IEEE754, @IEEE754]. The fact that almost all real numbers cannot be represented exactly in floating-point precision leads to important consequences for the numerical stability of any algorithm and the growth of numerical round-off error.

To our knowledge, we present the first publicly available Wisdom-Holman integrator that is unbiased, i.e. the errors are random and uncorrelated. This leads to a very slow error growth. For sufficiently small timesteps, we achieve Brouwer’s law, i.e., the energy error grows as 