Compute probability of a particular ordering of Multivariate normal distribution

Question

There are $m$ normally distributed, independent random variables $N_0,...,N_m$ with distinct means $\mu_1, \ldots \mu_m$ and standard deviations $\sigma_1, \ldots, \sigma_m$. Then, we get a permutation of the numbers $\{1, \ldots, m\}$. How can we efficiently compute, numerically, the (log) probability of observing the random variables in same ordering as this permutation?

An example:

1. we have four independent random variables $N1,N2,N3,N4$, all with different means and variances.
2. We are given the permutation $(3, 1, 2, 4)$.
3. What's $Pr(N3>N1>N2>N4)$?

I can get the answer by calculating the following formula.

$$\int_{-\infty}^\infty \int_{n_4}^\infty \int_{n_2}^\infty \int_{n_1}^\infty p(n_1)p(n_2)p(n_3)p(n_4)\ dn_3 dn_1 dn_2 dn_4$$

However, it is very difficult to analytically calculate this formula.

How do you think that it is better to calculate such multivariate normal CDF with higher accuracy?

I feel that having variables in this integration range is very troublesome.

I welcome any hints on calculation methods and means.

I do not know the details of the quasi-Monte Carlo method,,, Is it possible to apply quasi-Monte Carlo integration even when including variables in the integration interval like this?