Implementing PostGIS query using sf in R to make ST_Intersection faster

Question

I'm trying to implement this query that speeds up st_intersections by using faster functions, namely ST_Contains and ST_Within, to do ST_Intersection only where absolutely necessary, from PostGIS to R

SELECT CASE 
   WHEN ST_Contains(a.geom, b.geom)
     THEN b.geom 
   WHEN ST_Within(a.geom, b.geom)
     THEN a.geom 
   ELSE ST_Intersection(a.geom, b.geom)
   END ...
FROM a JOIN b
ON ST_Intersects(a.geom, b.geom)

what I have so far:

library(sf)
library(dplyr)
library(raster)

# anp_fed = raster::getData(country = "MEX", level = 1) %>% st_as_sf() %>% filter(NAME_1 == "Veracruz")
# veg     = read_sf("usv250ks4gw.shp") 
     # downloadable from http://www.conabio.gob.mx/informacion/gis/maps/geo/usv250ks4gw.zip

contienen  <- st_contains(anp_fed, veg, sparse = T)
contenidas <- st_within(anp_fed, veg, sparse = T)

contienen_list <- c()
for(i in 1:dim(contienen)[1]) { 
  contienen_list <- c(contienen[[i]], contienen_list)
}
contienen_list <- sort(contienen_list)

contenidas_list <- c()
t_contenidas <- t(contenidas)
for(j in 1:dim(t_contenidas)[1]) { 
  contenidas_list <- c(t_contenidas[[j]], contenidas_list)
}
contenidas_list <- sort(contenidas_list)


all_in_one <- rbind(st_join(anp_fed[contenidas_list,], veg),
      st_join(veg[contienen_list,], anp_fed) %>% as.data.frame() %>%
        dplyr::select(names(anp_fed)[-length(names(anp_fed))], names(veg)) %>% st_as_sf(),
      st_intersection(anp_fed[ifelse(length(contenidas_list) == 0, 1:dim(anp_fed)[1], -contenidas_list),],
                      veg[-contienen_list,])
      )

# system.time(source("postgis_to_sf.R"))     =  491 s                     
# system.time(st_intersection(anp_fed, veg)) =  666 s

I'm trying to do this without calling rpostgis because I want to do it directly on files in my system without having to put them into Postgresql; I think there's a lot of space for improvement, maybe there's a more straightforward way of translating sgbp to a vector to filter sf objects, and also the re-ordering of columns of the second st_join() call, transforming to data.frame; I'm trying to make a more concise code; note also the evaluation for empty sf objects needed (it should be done in every object, if it is to become a function)

Execution with st_intersection() takes 666 s, while the optimized version takes 491 s, almost 3 mins less.