Take a uniform on the surface of a d-dimensional sphere. Then multiply it by sqrt(chi2(df=d)). You'll get a MVN(0,I).

[swihart]

library(uniformly)

a1 <- -2.4
a2 <- 0.6
b1 <- -4
b2 <-  4
c1 <- -3
c2 <-  3

no.sims <- 1e6
sims <- runif_on_sphere(no.sims, d=3)
rchi <- sqrt(rchisq(no.sims,df=3))
prodd <- sims*rchi
mean(prodd[,1]>a1 & prodd[,1]<a2 & prodd[,2] <b2 & prodd[,2]>b1 & prodd[,3] < c2 & prodd[,3]>c1)
mvtnorm::pmvnorm(lower=c(a1,b1,c1),upper=c(a2,b2,c2))[1]


a1 <- -0.4
a2 <-  0.6
b1 <- -1.96
b2 <-  1.96

no.sims <- 3e6
sims <- runif_on_sphere(no.sims, d=2)
rchi <- sqrt(rchisq(no.sims,df=2))
prodd <- sims*rchi
mean(prodd[,1]>a1 & prodd[,1]<a2 & prodd[,2] <b2 & prodd[,2]>b1)
mvtnorm::pmvnorm(lower=c(a1,b1),upper=c(a2,b2))[1]


### THIS! for univariate ###
a1 <- 0.4
a2 <- 0.6

no.sims <- 3e6
sims <- runif_on_sphere(no.sims, d=1); table(sims)
rchi <- sqrt(rchisq(no.sims,df=1))
prodd <- sims*rchi
mean(prodd>a1 & prodd<a2)
pnorm(a2) - pnorm(a1)


### NOT-this! for univariate ###
a1 <- 0.4
a2 <- 0.6

no.sims <- 1e3
sims <- runif(no.sims, 0,1)
rchi <- sqrt(rchisq(no.sims,df=1))
prodd <- sims*rchi
mean(prodd>a1 & prodd<a2)
pnorm(a2) - pnorm(a1)