Multivariate normal cumulative distribution

state-of-the-art algorithm for computing the multivariate normal cdf in high dimensions

Computes the probability Pr(l<X<u), where 'X' is a zero-mean multivariate normal vector with covariance 'Sig'.
In high dimensions, this algorithm is vastly superior to the one in Matlab's statistics toolbox, see example.
Reference: Z. I. Botev (2015), "The Normal Law Under Linear Restrictions: Simulation and Estimation via Minimax Tilting", submitted to JRSS(B)