is in principle applicable to any precision. However, we work exclusively in double floating-point precision (64 bit) which is used on almost all modern CPUs.

Jacobi Coordinate Transformations {#sec:jacobi}
---------------------------------

The evolution under the effect of the interaction Hamiltonian is most efficiently done in Cartesian coordinates. On the other hand, the evolution of the Kepler Hamiltonian is easier in Jacobi coordinates. We thus need an efficient way to convert to and from Jacobi coordinates.

Luckily, the conversion from Cartesian to Jacobi coordinates and back can be done efficiently in $\mathcal{O}(N)$. We construct the algorithms from the definitions above and list them here in pseudo code. As before, primes denote Jacobi coordinates, Note that these algorithms work even if some of the bodies are test particles with $m_i=0$ (for $i\neq 0$). To convert from Cartesian to Jacobi coordinates:

$\mathbf{R} \gets m_0 \cdot \mathbf{r}_0$ $\mathbf{r}'_i \gets \mathbf{r}_i - \mathbf{R}/M_{i-1}$ $\mathbf{R}  \gets \mathbf{R}  \cdot (1 +m_i/M_{i-1}) + m_i\cdot\mathbf{r}'_i$ $\mathbf{r}'_0 \gets \mathbf{R} / M_{N-1}$

Similarly, we construct the algorithm to convert back from Jacobi to Cartesian coordinates as follows:

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

