Image simulation for scanning transmission electron microscopy at atomic resolution for samples with realistic dimensions can require very large computation times using existing simulation algorithms. We present a new algorithm named PRISM that combines features of the two most commonly used algorithms, namely the Bloch wave and multislice methods. PRISM uses a Fourier interpolation factor f that has typical values of 4–20 for atomic resolution simulations. We show that in many cases PRISM can provide a speedup that scales with f^4 compared to multislice simulations, with a negligible loss of accuracy. We demonstrate the usefulness of this method with large-scale scanning transmission electron microscopy image simulations of a crystalline nanoparticle on an amorphous carbon substrate. We also present a software package called Prismatic for parallelized simulation of image formation in STEM using both the PRISM and multislice methods. By distributing the workload between multiple CUDA-enabled GPUs and multicore processors, accelerations as high as 1000x for PRISM and 15x for multislice are achieved relative to traditional multislice implementations using a single 4-GPU machine.Prismatic is freely available both as an open-source CUDA/C++ package with a graphical user interface for Windows and OSX, and as a Python package.