From: Petr Savicky <savicky_at_cs.cas.cz>

Date: Mon, 15 Oct 2007 12:49:22 +0200

R-devel_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-devel Received on Mon 15 Oct 2007 - 10:57:08 GMT

Date: Mon, 15 Oct 2007 12:49:22 +0200

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.

R-devel_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-devel Received on Mon 15 Oct 2007 - 10:57:08 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Thu 25 Oct 2007 - 11:37:11 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel.
Please read the posting
guide before posting to the list.
*