`s2_union_agg()` is slow

The strategy of accumulating a union is correct but slow! I assume there is a faster way to go about this that I haven't found in the documentation yet.

library(s2)
library(sf)
#> Warning: package 'sf' was built under R version 4.0.5
#> Linking to GEOS 3.8.1, GDAL 3.2.0, PROJ 7.2.0

countries <- s2_data_countries()
countries_sf <- st_as_sf(countries)

bench::mark(
  s2_coverage_union_agg(countries),
  s2_union_agg(countries),
  st_union(countries_sf),
  check = F
)
#> # A tibble: 3 x 6
#>   expression                            min   median `itr/sec` mem_alloc
#>   <bch:expr>                       <bch:tm> <bch:tm>     <dbl> <bch:byt>
#> 1 s2_coverage_union_agg(countries) 100.89ms  101.7ms     9.59     61.7KB
#> 2 s2_union_agg(countries)             2.82s    2.82s     0.355    43.7KB
#> 3 st_union(countries_sf)           141.69ms 144.13ms     6.90    666.1KB
#> # … with 1 more variable: gc/sec <dbl>

^{Created on 2021-04-28 by the reprex package (v0.3.0)}

Start by reproducing the benchmark for s2_union_agg(countries) and compare it with s2_coverage_union_agg() and st_union(). Trace the s2_union_agg entry point and its underlying implementation; done means the union remains correct and the reported performance gap is substantially reduced, with benchmark coverage added or updated.