Reliable extraction of cosmological information from clustering measurements of galaxy surveys requires estimation of the error covariance matrices of observables. The accuracy of covariance matrices is limited by our ability to generate sufficiently large number of independent mock catalogs that can describe the physics of galaxy clustering across a wide range of scales. Furthermore, galaxy mock catalogs are required to study systematics in galaxy surveys and to test analysis tools. In this investigation, we present a fast and accurate approach for generation of mock catalogs for the upcoming galaxy surveys. Our method relies on low-resolution approximate gravity solvers to simulate the large scale dark matter field, which we then populate with halos according to a flexible nonlinear and stochastic bias model. In particular, we extend the textsc{patchy} code with an efficient particle mesh algorithm to simulate the dark matter field (the textsc{FastPM} code), and with a robust MCMC method relying on the textsc{emcee} code for constraining the parameters of the bias model. Using the halos in the BigMultiDark high-resolution $N$-body simulation as a reference catalog, we demonstrate that our technique can model the bivariate probability distribution function (counts-in-cells), power spectrum, and bispectrum of halos in the reference catalog. Specifically, we show that the new ingredients permit us to reach percentage accuracy in the power spectrum up to $ksim 0.4; ,h,{rm Mpc}^{-1}$ (within 5% up to $ksim 0.6; ,h,{rm Mpc}^{-1}$) with accurate bispectra improving previous results based on Lagrangian perturbation theory.