simp.d.dim <- function(n,d){
  y <- seq(0,5,5/(n-1))
  wsimp <- c(1,rep(c(4,2),(n-3)/2),4,1)
  simp <- ((5/(3*(n-1)))*sum((1/(1+y))*wsimp))^d
  return(simp)
}

simp.d.dim(5,1)


mc.d.dim <- function(n,d){
  u <- matrix(5*runif(n*d),ncol=d)
  mc.int <- (5^d)*mean(apply(1/(1+u),1,prod))
  mc.sd <- (5^d)*sd(apply(1/(1+u),1,prod))
  return(c(mc.int,mc.sd))
}

mc.d.dim(5^2,1)