Replacing the Random Number Generator in Stand Alone Library

Okay, so I am guessing everyone had the same response I initially did when
hearing that this RNG might not be so hot*...  as an alternate question, if
I submitted a patch to replace the current RNG with the twister, would it be
accepted?

RNGversion("1.6.0")

s <- replicate(1e6, max(rnorm(100)))
plot(density(s))
m <- matrix(runif(8e7),2)
plot(m[1,],m[2,], xlim=c(0,1e-3), pch=".")

-N

-----Original Message-----
From: Nigel Delaney [mailto:nigel.delaney at outlook.com] 
Sent: Thursday, October 10, 2013 5:04 PM
To: r-devel at r-project.org
Subject: Replacing the Random Number Generator in Stand Alone Library

Hi R-Developers,

I had a question about the random number generator used in the R StandAlone
Math Library.  The stand-alone library depends on the unif_rand() function
for most simulated values, and this function is provided in the sunif.c file
in the relevant directory.  At present, this program implements the
"Marsaglia-Multicarry" algorithm, which is described throughout the R
documentation as:

"A multiply-with-carry RNG is used, as recommended by George Marsaglia in
his post to the mailing list 'sci.stat.math'. It has a period of more than
2^60 and has passed all tests (according to Marsaglia). The seed is two
integers (all values allowed)."

However, I do not think this RNG actually passes all tests.   For example,
the Handbook of Computational Econometrics (illegal web copy at link below),
shows that it fails the mtuple test and gives an explicit example where it
leads to problems because it failed this test.  The mtuple test was
introduced by Marsaglia in 1985, and I gather he wrote his mailing list
comment that it "passes all tests" sometime after this, so I am not sure
what explains this distinction (though I am not sure if the mtuple test is
included in the diehard tests, which he may have been what he was referring
to).  However, there are clearly some areas where this PRNG runs in to
trouble (although the books example is better, another problem is that it
can't seem to simulate a value above (1/2)^1+(1/4)^4) after simulating a
value below 1e-6.

Given that the Mersenne Twister seems to be the standard for simulation
these days (and used as the default in R), it seems like it might be useful
to change the stand alone library so it also uses this routine.  I gather
this would be pretty easy to do by pulling this function from the RNG.c file
and moving it into the sunif.c file, and have a prototype of this.

However, I thought I would ask, is there a reason this hasn't been done?  Or
is it just a historical carry-over (pun intended I suppose).

Warm wishes,
Nigel

********************
Research Fellow
Massachusetts General Hospital / Broad Institute


Book link:
http://thenigerianprofessionalaccountant.files.wordpress.com/2013/04/handboo
k-of-computational-econometrics-belsley.pdf