  A boundary integral equation method for the 3-D Helmholtz equation in
multilayered media with many quasi-periodic layers is presented. Compared with
conventional quasi-periodic Green's function method, the new method is robust
at all scattering parameters. A periodizing scheme is used to decompose the
solution into near- and far-field contributions. The near-field contribution
uses the free-space Green's function in an integral equation on the interface
in the unit cell and its immediate eight neighbors; the far-field contribution
uses proxy point sources that enclose the unit cell. A specialized high-order
quadrature is developed to discretize the underlying surface integral operators
to keep the number of unknowns per layer small. We achieve overall linear
computational complexity in the number of layers by reducing the linear system
into block tridiagonal form and then solving the system directly via block LU
decomposition. The new solver is capable of handling a 100-interface structure
with 961.3k unknowns to $10^{-5}$ accuracy in less than 2 hours on a desktop
workstation.
