B <- function(z,n) sqrt(1+n)*exp(-(z^2/2)/(1+1/n))
p0 <- function(z,n) (1 + B(z,n)^(-1))^(-1)
p <- c(0.1,0.05,0.01,0.001)
k <- c(1,5,10,20,50,100,1000)
nr <- length(p)
nc <- length(k)
z <- rep(0,nr)
cat("P-value  n\\z")
  for (j in 1:nc) 
    cat("  ",formatC(k[j],width=4,digits=0,format="f"))
cat("\n")
for (i in 1:4){
  z[i] <- qnorm(1-p[i]/2)
  cat(formatC(p[i],width=4,digits=3,format="f")," ",
      formatC(z[i],width=4,digits=3,format="f"))
  for (j in 1:nc) 
    cat(" ",formatC(p0(z[i],k[j]),width=4,digits=3,format="f"))
  cat("\n")
}
cat("\n")
cat("P-value   z    Bound on B   Bound on p0\n")
for (i in 1:4){
  z[i] <- qnorm(1-p[i]/2)
  cat(formatC(p[i],width=4,digits=3,format="f")," ",
      formatC(z[i],width=4,digits=3,format="f"))
  cat("   ",formatC(exp(-z[i]^2/2),width=4,digits=3,format="f"))
  cat("       ",
      formatC((1+exp(z[i]^2/2))^(-1),width=4,digits=3,format="f"))
  cat("\n")
}