A computationally efficient algorithm is presented for exact simulation of the stochastic time evolution of spatially homogeneous aggregation-fragmentation processes featuring multiple components or conservation laws. The algorithm can predict the average size and composition distributions of aggregating particles as well as their fluctuations, regardless of the functional form (e.g., composition dependence) of the aggregation or fragmentation kernels. Furthermore, it accurately predicts the complete time evolutions of all moments of the size and composition distributions, even for systems that exhibit gel transitions. We demonstrate the robustness and utility of the algorithm in case studies of linear and branched polymerization processes, the last of which is a two-component process. These simulation results provide the stochastic description of these processes and give new insights into their gel transitions, fluctuations, and long-time behavior when deterministic approaches to aggregation kinetics may not be reliable.