  Context. Due to our increasing knowledge on the Galactic and stellar
neighborhood of the Solar System, modern long-period comet motion studies have
to take into account both stellar perturbations and the overall Galactic
potential. Aims. Our aim is to propose algorithms and methods to perform
numerical integration of a Solar System small body equations of motion much
faster and at the same time with greater precision. Methods. We propose a new
formulation of the equations of motion formulated in the Solar System
barycentric frame but accurately accounting for the differential perturbations
caused by the Galactic potential. To use these equations effectively we provide
numerical ephemerides of the Galactic positions of the Sun and a set of
potential stellar perturbers. Results. The proposed methods offer the precision
higher by several orders of magnitude and simultaneously greatly reduce the
necessary CPU time. The application of this approach is presented with the
example of a detailed dynamical study of the past motion of comet C/2015 XY1.
