  We present a new boundary integral formulation for time-harmonic wave
diffraction from two-dimensional structures with many layers of arbitrary
periodic shape, such as multilayer dielectric gratings in TM polarization. Our
scheme is robust at all scattering parameters, unlike the conventional
quasi-periodic Green's function method which fails whenever any of the layers
approaches a Wood anomaly. We achieve this by a decomposition into near- and
far-field contributions. The former uses the free-space Green's function in a
second-kind integral equation on one period of the material interfaces and
their immediate left and right neighbors; the latter uses proxy point sources
and small least-squares solves (Schur complements) to represent the remaining
contribution from distant copies. By using high-order discretization on
interfaces (including those with corners), the number of unknowns per layer is
kept small. We achieve overall linear complexity in the number of layers, by
direct solution of the resulting block tridiagonal system. For device
characterization we present an efficient method to sweep over multiple incident
angles, and show a $25\times$ speedup over solving each angle independently. We
solve the scattering from a 1000-layer structure with $3\times 10^5$ unknowns
to 9-digit accuracy in 2.5 minutes on a desktop workstation.
