I used Haskell to do this. Code is here I did not use standard `System.Random` library to generate random number but `mwc-random` package. I can’t say if speed is better or worse. I wanted to see how fast we converge to the value of `pi`. On x-axis, we have number of samples (N), and for each N, we compute the value of pi 50 times. This way we can talk about mean and variance in computed pi (which we plot on y-axis). It seems that for N >= 1000, we get a reasonable approximation of pi with low variance. For a 2-D sphere (aka circle: \$latex x_1^2 + x_2 ^ 2 = 1 \$ ), The convergence is following. For a 3-D sphere, (that is \$latex x_1 ^2 + x_2 ^2 + x_3 ^3 = 1\$), I get the following: Now I am bit more ambitious, I use a 10-D sphere: As you can see, we converge to a a reasonable value of pi with low variance but convergence is not improving monotonically as in previous cases. The mean starts from a very low value and get to a respectable value only after N > 1000 (in low-dimensions cases we were swinging around 3.1415). The answer: we need to sample minimum number of points from hyper-cube to get some reasonable estimate. How variance decreases with N (when N is significantly large, 100 repeats of each computation)? Following plots starts from 2D and goes onto 10 dimensions.           