  We present a solver for plane wave scattering from a periodic dielectric
grating with a large number $M$ of inclusions lying in each period of its
middle layer.Such composite material geometries have a growing role in modern
photonic devices and solar cells. The high-order scheme is based on boundary
integral equations, and achieves many digits of accuracy with ease. The usual
way to periodize the integral equation---via the quasi-periodic Green's
function---fails at Wood's anomalies. We instead use the free-space Green's
kernel for the near field, add auxiliary basis functions for the far field, and
enforce periodicity in an expanded linear system; this is robust for all
parameters. Inverting the periodic and layer unknowns, we are left with a
square linear system involving only the inclusion scattering coefficients.
Preconditioning by the single-inclusion scattering matrix, this is solved
iteratively in $O(M)$ time using a fast matrix-vector product. Numerical
experiments show that a diffraction grating containing $M=1000$ inclusions per
period can be solved to 9-digit accuracy in under 5 minutes on a laptop.
