Here is my function for a two-sample pooled binomial z-test based on NIST Engineering Statistics Handbook: (all the p’s are actually p^hat’s)

bin_test <- function(x1, n1, x2, n2) {
  p1 <- x1 / n1
  p2 <- x2 / n2
  p <- (x1 + x2) / (n1 + n2)
  z <- (p2 - p1) / sqrt(p*(1-p)*(1/n1 + 1/n2))
  z
}

Does this exist as an R built-in in stats or package? I could only find prop.test() which for 2-sample gives the chi-squared statistic which is almost the z-score squared.