We demonstrate optimization of optical metasurfaces over 105-106 degrees of freedom in two and three dimensions, 100-1000+ wavelengths (λ) in diameter, with 100+ parameters per λ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 λ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.