We describe a fully discrete high-order algorithm for simulating multiple scattering of electromagnetic waves in three dimensions by an ensemble of perfectly conducting scattering objects. A key component of our surface integral algorithm is high-order tangential approximation of the surface current on each obstacle in the ensemble. The high-order nature of the algorithm leads to relatively small numbers of unknowns, which allows us to use either a direct method or an iterative boundary decomposition method for simulations of multiple scattering. We demonstrate the algorithm using both of these techniques for near and well separated obstacles. Using a small computing cluster (with 20 processors), we simulate multiple scattering by up to 125 objects for frequencies in the resonance region, and by paired obstacles of diameter 20 to 30 times the incident wavelength. Many physically important problems, such as scattering by atmospheric aerosols or ice crystals, involve multiple scattering by ensembles of particles, each particle having its own distinct shape, but with all particles fitting a stochastic description with a small number of fixed parameters in local spherical coordinates. We demonstrate our algorithm for multiple scattering by ensembles of such unique particles, whose stochastic description corresponds to computer models of ice crystals and dust particles.