High-accuracy RNG I need to simulate a sample of size 100x3 of white noise (normal with null mean and unitary variance)
I'm actually using Mersenne-Twister but discussing with my Professor seems is not "very stable" for the size we need.
Since I'm going to perform it via a programming language, do you have any suggestion for other RNG? I also found this.
Thanks for the help :)
UPDATE
With "not very stable" I mean how pretty different the values are each run. rnorm should give me pretty independent elements with are not since covariance is not equal to Identity Matrix. I would need to have a sample that result pretty close to identity matrix, even if rounded. So, is there any way to discard samples until it's very close to Identity Matrix? What condition should I use?

 A: The fact that the variances are not exactly 1 and the covariances not exactly 0 is totally expected, especially with these sample sizes and does by no means indicate that the RNG is not working properly. (Hint: Try running cov(matrix(rnorm(1e6*3), 1e6))). Remember that the values specified in the RNG are population parameters of the distribution and not the empirical values of samples drawn from this population.
Although it is unclear from the question why such exact numbers are needed, they can be generated using the function mvrnorm from the MASS package specifying the argument empirical = TRUE. In contrast to rnorm which is used to sample from univariate normal distributions, mvrnorm is used to generate samples from multivariate normal distributions. Here is an example:
library(MASS)

set.seed(142857)
mu_vec <- c(0, 0, 0)  # The means
cov_mat <- diag(1, 3) # The covariance matrix
n <- 100              

x <- mvrnorm(n = n, mu = mu_vec, Sigma = cov_mat, empirical = TRUE)

cov(x)
              [,1]          [,2]          [,3]
[1,]  1.000000e+00 -2.250771e-16 -2.604286e-16
[2,] -2.250771e-16  1.000000e+00  6.497218e-17
[3,] -2.604286e-16  6.497218e-17  1.000000e+00

The variances and covariances of the simulated samples are now exactly 1 and 0, respectively, as desired.
A: All right, so it looks to me that your problem has nothing to do with the random number generator but rather with choosing a distribution that does not generate a lot of outliers (low kurtosis)
If there are no constraints preventing us from doing so, I would go for the uniform distribution:
matrix(runif(300, min=0, max=1), 100)

ADDED AFTER QUESTION EDIT: Since you can't change your distribution. I think your best chance is to save the output of first execution (setting variance to $x$) and add some random noise (normal with variance $1-x$ This way you get the desired mean and variance while getting some stability between executions. Choose the value of $x$ that better fits your needs
