N <- sum(n)
Z <- (n - N*p0) / sqrt(N * p0 * (1 - p0))

k <- length(n)
maxz <- max(abs(Z))
pvalue <- 2 * k * (1 - pnorm(maxz))