  We demonstrate optimization of optical metasurfaces over $10^5$--$10^6$
degrees of freedom in two and three dimensions, 100--1000+ wavelengths
($\lambda$) in diameter, with 100+ parameters per $\lambda^2$. In particular,
we show how topology optimization, with one degree of freedom per
high-resolution "pixel," can be extended to large areas with the help of a
locally periodic approximation that was previously only used for a few
parameters per $\lambda^2$. In this way, we can computationally discover
completely unexpected metasurface designs for challenging multi-frequency,
multi-angle problems, including designs for fully coupled multi-layer
structures with arbitrary per-layer patterns. Unlike typical metasurface
designs based on subwavelength unit cells, our approach can discover both sub-
and supra-wavelength patterns and can obtain both the near and far fields.
