  The general consensus in the N-body community is that statistical results of
an ensemble of collisional N-body simulations are accurate, even though
individual simulations are not. A way to test this hypothesis is to make a
direct comparison of an ensemble of solutions obtained by conventional methods
with an ensemble of true solutions. In order to make this possible, we wrote an
N-body code called Brutus, that uses arbitrary-precision arithmetic. In
combination with the Bulirsch--Stoer method, Brutus is able to obtain converged
solutions, which are true up to a specified number of digits. We perform
simulations of democratic 3-body systems, where after a sequence of resonances
and ejections, a final configuration is reached consisting of a permanent
binary and an escaping star. We do this with conventional double-precision
methods, and with Brutus; both have the same set of initial conditions and
initial realisations. The ensemble of solutions from the conventional
simulations is compared directly to that of the converged simulations, both as
an ensemble and on an individual basis to determine the distribution of the
errors. We find that on average at least half of the conventional simulations
diverge from the converged solution, such that the two solutions are
microscopically incomparable. For the solutions which have not diverged
significantly, we observe that if the integrator has a bias in energy and
angular momentum, this propagates to a bias in the statistical properties of
the binaries. In the case when the conventional solution has diverged onto an
entirely different trajectory in phase-space, we find that the errors are
centred around zero and symmetric; the error due to divergence is unbiased, as
long as the time-step parameter, eta <= 2^(-5) and when simulations which
violate energy conservation by more than 10% are excluded.
