We thoroughly tested the conversions to and from Jacobi coordinates to ensure they are unbiased. This task turns out to be much harder than we na[ï]{}vely expected. As an example, consider the following algorithm which is formally equivalent to the above but numerically much less stable.

$\mathbf{R} \gets 0$ $\mathbf{r}_i \gets \mathbf{r}'_0 + M_{i-1}/M_i \cdot \mathbf{r}'_i -\mathbf{R}$ $\mathbf{R}  \gets \mathbf{R}  + m_i/M_i \cdot \mathbf{r}'_i$ $\mathbf{r}_0 \gets \mathbf{r}'_0 - \mathbf{R}$

In the above algorithm, we access $\mathbf{r}_0'$ multiple times and have to do a subtraction in the last step. This significantly promotes error propagation and leads to floating-point errors that can be orders of magnitudes higher than in the other implementation. After many timesteps, this leads to a linear secular growth in the energy error.

Implementation of Newton’s Method {#sec:newton}
---------------------------------

To solve Kepler’s equation (Eq. \[eq:kepler\]) for $X$, we need to use an iterative scheme. We now describe our implementation of Newton’s method in floating-point arithmetic. The straightforward implementation is an iteration loop that terminates when the change to $X$ is small, e.g.,

$X \gets initial\; guess$ $dX \gets -f(X)/f'(X)$ $X \gets X + dX$ .

Here, ${\ensuremath{\mathit{eps}}}$ is a small number just above machine precision, typically ${\ensuremath{\mathit{eps}}}\sim10^{-15}$. We use 