Monte-Carlo simulations of the fracture of elastic unidirectional model fibre composites are an important tool to understand composite reliability. On account of being computationally intensive, fracture simulations reported in the literature have been limited to simulation patches comprised of a few thousand fibres. While these limited patch sizes suffice to capture the dominant failure event when the fibre strength variability is low (synthetic fibres), they suffer from edge effects when the fibre strength variability is high (natural fibres). On the basis of recent algorithmic developments based on Fourier acceleration, a novel bisection based Monte Carlo failure simulation algorithm is presently proposed. This algorithm is used to obtain empirical strength distributions for model composites comprised of up to 2 20≈ 10 6 fibres, and spanning a wide range of fibre strength variabilities. These simulations yield empirical weakest-link strength distributions well into the lower tail. A stochastic model is proposed for the weakest-link event. The strength distribution predicted by this model fits the empirical distributions for any fibre strength variability. © 2019, Springer Nature B.V.