A method for generating theoretical breakage distribution functions for multiple particle breakage is presented. It starts with the joint probability function that accounts for all the child particles; it is then reduced to the marginal probability function commonly used in the breakage equation. This method is flexible enough to allow the user to choose the number of child particles and the functional form to be used. The method is demonstrated with both product and summation functions with a power-law form. To facilitate the use of these theoretical functions for statistical analyses, a companion discretized breakage equation is developed. The new equation guarantees the conservation of mass and correct prediction of the total number of particles despite discretization. It is easy to use because it is a set of ordinary differential equations and applicable to both equal-size and geometric-size intervals. Simulation results show that different breakage distribution functions coupled with different breakage rates can produce almost indistinguishable particle-size distributions signifying the need for further work in this area.