  We introduce FastPM, a highly-scalable approximated particle mesh N-body
solver, which implements the particle mesh (PM) scheme enforcing correct linear
displacement (1LPT) evolution via modified kick and drift factors. Employing a
2-dimensional domain decomposing scheme, FastPM scales extremely well with a
very large number of CPUs. In contrast to COmoving-LAgrangian (COLA) approach,
we do not require to split the force or track separately the 2LPT solution,
reducing the code complexity and memory requirements. We compare FastPM with
different number of steps ($N_s$) and force resolution factor ($B$) against 3
benchmarks: halo mass function from Friends of Friends halo finder, halo and
dark matter power spectrum, and cross correlation coefficient (or
stochasticity), relative to a high resolution TreePM simulation. We show that
the modified time stepping scheme reduces the halo stochasticity when compared
to COLA with the same number of steps and force resolution. While increasing
$N_s$ and $B$ improves the transfer function and cross correlation coefficient,
for many applications FastPM achieves sufficient accuracy at low $N_s$ and $B$.
For example, $N_s=10$ and $B=2$ simulation provides a substantial saving (a
factor of 10) of computing time relative to $N_s=40$, $B=3$ simulation, yet the
halo benchmarks are very similar at $z=0$. We find that for abundance matched
halos the stochasticity remains low even for $N_s=5$. FastPM compares well
against less expensive schemes, being only 7 (4) times more expensive than 2LPT
initial condition generator for $N_s=10$ ($N_s=5$). Some of the applications
where FastPM can be useful are generating a large number of mocks, producing
non-linear statistics where one varies a large number of nuisance or
cosmological parameters, or serving as part of an initial conditions solver.
