A general recursion algorithm is described for calculating kinematical diffraction intensities from crystals containing coherent planar faults. The method exploits the self-similar stacking sequences that occur when layers stack non-deterministically. Recursion gives a set of simple relations between average interference terms from a statistical crystal, which can be solved as a set of simultaneous equations. The diffracted intensity for a polycrystalline sample is given by the incoherent sum of scattered intensities over an ensemble of crystallites. The relations between this and previous approaches, namely the Hendricks-Teller matrix formulation, the difference equation method, the summed series formula of Cowley, and Michalski's recurrence relations between average phase factors, are discussed. Although formally identical to these previous methods, the present recursive description has an intuitive appeal and proves easier to apply to complex crystal structure types. The method is valid for all types of planar faults, can accommodate long-range stacking correlations, and is applicable to crystals that contain only a finite number of layers. A FORTRAN program DIFFaX, based on this recursion algorithm, has been written and used to simulate powder X-ray (and neutron) diffraction patterns and single crystal electron (kinematical) diffraction patterns. Calculations for diamond-lonsdaleite and for several synthetic zeolite systems that contain high densities of stacking faults are presented as examples.