Swarms of Stokeslets have previously been shown to be effective for simulating three-dimensional, free-surface, buoyancy-driven drop flows at unit viscosity ratio, including interfacial rupture and pass-through phenomena. This article presents an efficient marker-and-cell/fast Fourier transform (MAC-FFT) algorithm that yields N log N scaling of the operations, thereby enabling accurate simulations involving millions of particles. The formerly separate steps of regularization of the Green's function (commonly used in vortex blob methods for inviscid flow) and discretization are unified by regularizing the Stokeslet over the cubical cells that form the underlying gridas opposed to the spherical blobs used in the past. A piecewise-constant drop phase function, obtained by simple binning of the particles, thereby yields second-order quadrature of the Green's function to obtain the velocity field. An iterative cascade algorithm (adapted from Daubechies wavelets) allows the Stokes ocubeleto field to be calculated very efficiently. A similar cubelet approach is used for the cohesive forces that mimic interfacial tension. A user-friendly library of Fortran 95 subroutines (DropLib, www.droplib.org) has been developed to carry out the simulations and visualize drop shape evolutions in three dimensions.