  We present a novel hierarchical formulation of the fourth-order forward
symplectic integrator and its numerical implementation in the GPU-accelerated
direct-summation N-body code FROST. The new integrator is especially suitable
for simulations with a large dynamical range due to its hierarchical nature.
The strictly positive integrator sub-steps in a fourth-order symplectic
integrator are made possible by computing an additional gradient term in
addition to the Newtonian accelerations. All force calculations and kick
operations are synchronous so the integration algorithm is manifestly
momentum-conserving. We also employ a time-step symmetrisation procedure to
approximately restore the time-reversibility with adaptive individual
time-steps. We demonstrate in a series of binary, few-body and million-body
simulations that FROST conserves energy to a level of $|\Delta E / E| \sim
10^{-10}$ while errors in linear and angular momentum are practically
negligible. For typical star cluster simulations, we find that FROST scales
well up to $N_\mathrm{GPU}^\mathrm{max}\sim 4\times N/10^5$ GPUs, making direct
summation N-body simulations beyond $N=10^6$ particles possible on systems with
several hundred and more GPUs. Due to the nature of hierarchical integration
the inclusion of a Kepler solver or a regularised integrator with
post-Newtonian corrections for close encounters and binaries in the code is
straightforward.
