Monte Carlo simulations of the failure of unidirectional fibre composites typically require numerous evaluations of the stress-state in partially damaged composite patches. In a simulated composite patch comprised of N fibres, of which Nb fibres are broken in a common cross-sectional plane transverse to the fibre direction, the stress overloads in the intact fibres are given by the weighted superposition of the unit break solutions associated with each of the breaks. Determining the weights involves solving Nb linear equations, and determining overloads in the intact fibres requires matrix-vector multiplication. These operations require O(Nb3), and O(NNb) floating point operations, respectively. These costs become prohibitive for large N, and Nb; they limit Monte Carlo failure simulations to composite patches of only a few thousand fibres. In the present work, a fast algorithm to determine the overloads in a partially damaged composite, requiring O(Nb1/3NlogN) floating point operations, is proposed. This algorithm is based on the discrete Fourier transform. The efficiency of the proposed method derives from the computational simplicity of weighted superposition in Fourier space. Computations of the stress state ahead of large circular clusters of breaks in composite patches comprised of about one million fibres are used to demonstrate the efficiency of the proposed algorithm. © 2018, Springer Science+Business Media B.V., part of Springer Nature.