  Large-scale metasurfaces promise nanophotonic performance improvements to
macroscopic optics functionality, for applications from imaging to analog
computing. Yet the size scale mismatch of centimeter-scale chips versus
micron-scale wavelengths prohibits use of conventional full-wave simulation
techniques, and has necessitated dramatic approximations. Here, we show that
tailoring "fast direct" integral-equation simulation techniques to the form
factor of metasurfaces offers the possibility for accurate and efficient
full-wave, large-scale metasurface simulations. For cylindrical
(two-dimensional) metasurfaces, we demonstrate accurate simulations whose
solution time scales \emph{linearly} with the metasurface diameter. Moreover,
the solver stores compressed information about the simulation domain that is
reusable over many design iterations. We demonstrate the capabilities of our
solver through two designs: first, a high-efficiency, high-numerical-aperture
metalens that is 20,000 wavelengths in diameter. Second, a high-efficiency,
large-beam-width grating coupler. The latter corresponds to millimeter-scale
beam design at standard telecommunications wavelengths, while the former, at a
visible wavelength of 500 nm, corresponds to a design diameter of 1 cm, created
through full simulations of Maxwell's equations.
