[Rd] all zeroes in Mersenne Twister state may remain undetected

Petr Savicky savicky at cs.cas.cz
Mon Oct 15 12:49:22 CEST 2007


The function runif(n) contains a protection against a corrupted .Random.seed.
Besides other things, it tests whether at least one of the numbers 
.Random.seed[3:626] is not zero. If all(.Random.seed[3:626]==0), then
the internal Mersenne Twister state is regenerated using current time,
since a zero state is a fixed point of the recurrence and, hence, produces
a constant sequence.

However, the condition any(.Random.seed[3:626]!=0) does not imply that
the internal Mersenne Twister state is indeed not zero, since only the
most significant bit of .Random.seed[3] belongs to the internal state.
Hence, the number of bits in the state of Mersenne Twister
is 624*32 - 31 = 19937, which explains the period 2^19937-1.

For example, if .Random.seed[3] == 1, we always have the condition
  any(.Random.seed[3:626]!=0)
satisfied, but the internal state may still be effectively zero. An 
example of such a situation is
  RNGkind("default")
  .Random.seed[3:626] <- as.integer(0)
  .Random.seed[3] <- as.integer(1)
  x <- runif(10000)
  all(x==x[1]) # TRUE
  length(unique(x)) # 1
  all(.Random.seed[3:626]==0) # TRUE
Here, the internal state was effectively zero, but this fact was not
detected, since some (unimportant) bits of .Random.seed[3] were not zero.

On the contrary, if also .Random.seed[3]==0, then the internal state
is regenerated and the output of runif() becomes non constant:
  RNGkind("default")
  .Random.seed[3:626] <- as.integer(0)
  x <- runif(10000)
  all(x==x[1]) # FALSE
  length(unique(x)) # 10000
  all(.Random.seed[3:626]==0) # FALSE

The following patch to FixupSeeds corrects the detection of zero state:

--- R-devel_2007-10-14-orig/src/main/RNG.c	2007-09-02 07:49:35.000000000 +0200
+++ R-devel_2007-10-14-FixupSeeds/src/main/RNG.c	2007-10-15 04:33:52.988060624 +0200
@@ -181,7 +181,9 @@
 	 /* No action unless user has corrupted .Random.seed */
 	if(I1 <= 0) I1 = 624; 
 	/* check for all zeroes */
-	for (j = 1; j <= 624; j++)
+	notallzero = ((RNG_Table[RNG_kind].i_seed[1] & 0x80000000) != 0);
+	if(!notallzero)
+	for (j = 2; j <= 624; j++)
 	    if(RNG_Table[RNG_kind].i_seed[j] != 0) {
 		notallzero = 1;
 		break;

Petr Savicky.



More information about the R-devel mailing list