An algorithm is presented for rapid simulation of high-intensity focused ultrasound (HIFU) fields. Essentially, the method combines ray tracing with Monte Carlo integration to evaluate the Rayleigh-Sommerfeld integral. A large number of computational particles, phonons, are distributed among the elements of a phase-array transducer. The phonons are emitted into random directions and are propagated along trajectories computed with the ray tracing method. As the simulation progresses, an improving stochastic estimate of the acoustic field is obtained. The method can adapt to complicated geometries, and it is well suited to parallelization. The method is verified against reference simulations and pressure measurements from an ex vivo porcine thoracic tissue sample. Results are presented for acceleration with graphics processing units (GPUs). The method is expected to serve in applications, where flexibility and rapid computation time are crucial, in particular clinical HIFU treatment planning.