Theorem 5.3.1 in Larsen and Marx gives the large-sample Wald formula. The following function gives the score interval.
score <- function(n,x,level=.95) {
phat <- x/n
z <- -qnorm((1-level)/2)
moe <- ((z/sqrt(n))*sqrt(phat*(1-phat) + z^2/(4*n)))/(1+z^2/n)
center <- (phat+z^2/(2*n))/(1+z^2/n)
lcl <- center - moe
ucl <- center + moe
c(lcl,ucl)
}