Re: Issues with rand()

[Daniel J S. <dan...@ie...>]

Getting back to the 31-bit vs. 32-bit issue.  Here is the code of the algorithm:

  /* Generate pseudo random integers */
  k = seed1 / 53668L;
  seed1 = Xa1 * (seed1 - k * 53668L) - k * 12211;
  if (seed1 < 0)
	seed1 += Xm1;
  k = seed2 / 52774L;
  seed2 = Xa2 * (seed2 - k * 52774L) - k * 3791;
  if (seed2 < 0)
	seed2 += Xm2;
  z = seed1 - seed2;
  if (z < 1)
	z += (Xm1 - 1);

which has signed long values for the seeds.  So 017777777777L is correct.  However, I think the original comment means that negative values for the seeds should be allowed as well.  BUT, gnuplot is using negative values to mean reset the generator.  So we can't use negative values.  Furthermore, does it make sense that negative seeds should be allowed?

Well, let's verify that negative values for the seeds is sane and that the random number generator doesn't do something unexpected for negative values.  In the case of seed1, the most negative value results when seed1 is the largest value for which the division has no remainder, i.e., k = 40014 and seed1 = 2147471352.  In that case

seed1 = Xa1 * (seed1 - k * 53668L) - k * 12211;
     = Xa1 * (2147471352 - 40014 * 53668) - 40014 * 12211
     = Xa1 * 0 - 488610954
     = -488610954

then

if (seed1 < 0)
	seed1 += Xm1;

or

seed1 = -488610954 + 2147483563
     = 1658872609

In other words, in normal operation of the congruential random number generator, the seed1 value doesn't ever become negative after a pass of the algorithm.

Let's do the same check on seed2.  The most negative intermediate value happens for k = 40692 and seed2 = 2147479608, i.e.,

seed2 = Xa2 * (seed2 - k * 52774L) - k * 3791;
     = Xa2 * (2147479608 - 40692 * 52774) - 40692 * 3791
     = -154263372

if (seed2 < 0)
  seed2 += Xm2;

or

seed2 = -154263372 + 2147483399
     = 1993220027

So, again, if seed2 starts out positive, there is no way it can become negative after any iteration.

Not having the original paper in front of me, whether negative input seeds are allowed, I can't be sure.  But my feeling is that negative values should not be allowed input values because the normal operation has seed values which only land in the positive integer range.  The comment in the code is not a valid one, in my opinion.

If that is the case, then there is a bug in the code which hasn't been mentioned yet, e.g., rand({5,-20}).  That is, setting seed 1 greater than zero and seed 2 less than zero is not currently disallowed by the gnuplot code and could potentially put the generator in a strange state it normally would never see.

Ethan, please give the attached patch file a try.  I've included more detailed illegal seed values including rejection of any non-integer seeds and any inputs where seed2 is negative.  Also, the internal seed is not reset until after a valid reset input sequence is confirmed.  Also, I beefed up the documentation including a short sentence about the random number generation algorithm.

Here are some example outputs:

Terminal type set to 'x11'
gnuplot> print rand(2**31-1)
0.996865893971871
gnuplot> print rand(2**31)
         Illegal seed value

gnuplot> print rand(0.5)
         Illegal seed value

gnuplot> print rand(-1.5)
         Illegal seed value

gnuplot> print rand(-1)
0.222457440974512
gnuplot> print rand({-1,5})
         Illegal seed value

gnuplot> print rand({5,-1})
         Illegal seed value

gnuplot> print rand({5,0})
0.999998420858381
gnuplot> print rand({5,2147483648})
         Illegal seed value

gnuplot> print rand({5,0.5})
         Illegal seed value

I ran the above on the current CVS version of gnuplot and realized that 2^31-1 is currently not a valid seed.  I'm not sure why that should be disallowed.  In the attached patch, 2^31-1 is a valid seed value.

Dan