I find the previous code for the multivariate gaussian totally crap. https://xinxiwang.wordpress.com/2011/03/23/multivariate-gaussian/

Since the covariance matrix of gaussian distribution is symmetric positive definite which can be decomposed into LL’ by Cholesky decomposition, where L is a lower triangular matrix. Cholesky decomposition is much faster and stable than the general LU decomposition algorithm.

Another place can be improved is the determinant. After the LU decomposition is done, we can easily get the determinant by multiply together the diagonal elements of L and U, which can be done together with inversion in one function.