any size (at the price of increasing computational effort), it is sufficient to reduce it well below the level of the estimation error. We will now first assess the approximation error alone and then consider the total error.

Accuracy of the Approximation {#sec:accur}
-----------------------------

We estimate the accuracy of the approximated forces for various choices of the parameters controlling the algorithm and for two typical astrophysical situations: a spherical Hernquist [@Hern] model galaxy, which has density and force per unit mass $$\label{hernquist}
  \rho({{\mbox{\boldmath$x$}}}) = {{\sf M}_{\rm tot} \over 2\pi}
  {r_0\over |{{\mbox{\boldmath$x$}}}| \big(r_0+|{{\mbox{\boldmath$x$}}}|\big)^3},\qquad
  {{\mbox{\boldmath$F$}}}({{\mbox{\boldmath$x$}}})  = - {{{\mbox{\boldmath$x$}}}\over|{{\mbox{\boldmath$x$}}}|}
  { G{\sf M}_{\rm tot}\over\big(r_0+|{{\mbox{\boldmath$x$}}}|\big)^2},$$ and a group of five such galaxies at different positions and with various scale radii $r_0$. In either case, we sample $N=10^4$, $10^5$, and $10^6$ bodies. We use standard Plummer softening, where $g(r)=G/\sqrt{r^2+\epsilon^2}$, with softening lengths $\epsilon$ chosen such as to minimize the estimation error [@Db]. In order to single out the approximation error, we compare with the forces obtained by a computation via direct summation (in double precision; in case of $N=10^6$ for the first $10^5$ bodies only). As measure for the relative error, we compute for each body [@CM] $$\label{error}
  \varepsilon \equiv |a_{\rm approx} - a_{\rm direct}| / 