There are many ways to simulate a multivariate gaussian distribution assuming that you can simulate from independent univariate normal distributions.
One of the most popular method is based on the Cholesky decomposition.
Let’s see how Rcpp and Armadillo perform on this task.

The easiest way to perform a Cholesky distribution in R is to use the chol function in R which interface some fast LAPACK routines.

We will also use MASS:mvrnorm which implements it differently:

The RcppArmadillo function outperforms the MASS implementation and the naive R code, but more surprisinugly mvrnormR is slightly faster than mvrnorm in this benchmark.

To be fair, while digging into the MASS::mvrnorm code it appears that there are few code sanity checks ( such as the positive definiteness of Sigma ).