  Compared to other symplectic integrators (the Wisdom and Holman map and its
higher order generalizations) that also take advantage of the hierarchical
nature of the motion of the planets around the central star, our methods
require solving implicit equations at each time-step. We claim that, despite
this disadvantage, FCIRK16 is more efficient than explicit symplectic
integrators for high precision simulations thanks to: (i) its high order of
precision, (ii) its easy parallelization, and (iii) its efficient
mixed-precision implementation which reduces the effect of round-off errors. In
addition, unlike typical explicit symplectic integrators for near Keplerian
problems, FCIRK16 is able to integrate problems with arbitrary perturbations
(non necessarily split as a sum of integrable parts). We present a novel
analysis of the effect of close encounters in the leading term of the local
discretization errors of our integrator. Based on that analysis, a mechanism to
detect and refine integration steps that involve close encounters is
incorporated in our code. That mechanism allows FCIRK16 to accurately resolve
close encounters of arbitrary bodies. We illustrate our treatment of close
encounters with the application of FCIRK16 to a point mass Newtonian 15-body
model of the Solar System (with the Sun, the eight planets, Pluto, and five
main asteroids) and a 16-body model treating the Moon as a separate body. We
also present some numerical comparisons of FCIRK16 with a state-of-the-art high
order explicit symplectic scheme for 16-body model that demonstrate the
superiority of our integrator when very high precision is required.
