Re: Issues with rand()

[Ethan M. <merritt@u.washington.edu>]

On Tuesday, December 28, 2010, Juhász Péter wrote:
> Dear gnuplot developers,
> 
> a recent thread[1] in the newsgroup made me play with the rand()
> function, and I've found some issues with it:
> 
> 1) rand() silently accepts string arguments:
> gnuplot> print rand("foo")
> 0.222457440974512
> 
> While this is not particularly bad in itself, but it's not consistent
> with the behavior of other numerical functions (which reject string
> arguments), and it's undocumented. Its effect is the same as that of
> rand(0) - which have misled the user in the thread[1], because he
> thought that rand("time") sets the random seed to the current time.

I am not terribly concerned that passing a string to a numeric function
produces a garbage result, but it would be easy enough to do the same
thing in specfun.c that is already done in internal.c and standard.c
  internal.c:#define pop(x) pop_or_convert_from_string(x)
  standard.c:#define pop(x) pop_or_convert_from_string(x)

 > 2) It is possible to lock the PRNG into a state where the returned
> numbers are not random at all:
> gnuplot> print rand(0.5)
> 0.999999999450207
> gnuplot> print rand(0)
> 0.999999999450207
> gnuplot> print rand(0)
> 0.999999999450207
> gnuplot> print rand(0)
> 0.999999999450207
> gnuplot> print rand(0)
> 0.999999999450207

Here I agree that the documentation should state that seeds are
required to be non-zero integers, although it is then tricky to
explain that the way to set both seed values is
   rand( i + j*{0,1} )

It's fair to consider acceptance of a seed value (0 < seed < 1) as a
bug, since it will immediately be converted to an integer value 0
which is not a valid seed.  This might be an adequate fix:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
--- gnuplot/src/specfun.c       2010-10-21 22:28:24.000000000 -0700
+++ gnuplot-cvs/src/specfun.c   2010-12-28 13:23:20.000000000 -0800
@@ -1098,7 +1098,7 @@ ranf(struct value *init)
 
     /* Construct new seed values from input parameter */
     /* FIXME: Ideally we should allow all 64 bits of seed to be set */
-    if (real(init) > 0.0) {
+    if (real(init) > 1.0) {
        if (real(init) >= (double)(017777777777UL))
            int_error(NO_CARET,"Illegal seed value");
        if (imag(init) >= (double)(017777777777UL))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Note that the FIXME comment is out of date.  You can indeed set 64 bits
of seed by using rand(i + j*{0,1}) as noted above.


> 3) Consider the following plot:
> 
> set xrange [0:2**22]
> set samples 1000
> plot '+' u 1:(rand($1)) w l
> 
> The rand() function, when called with a nonzero argument, sets the seeds
> based on the argument and returns the first pseudo-random number from
> the sequence associated those seeds. As the plot shows, there is a
> rather obvious dependence between rand's argument and the returned
> number, in fact the dependence is linear if only the lower 21-or-so bits
> of the argument are considered.

Now I'm out of my depth.  
Pseudo-random number generation is fraught with pitfalls.

Is it necessarily a bad thing that the first number is related to the seed?  
The primary requirement is that successive numbers within a generated
sequence be sufficiently independent, or so I understand it.

A requirement that related seeds produce unrelated starting values
for the sequence seems like a quite different thing.  I can imagine that 
such a property may be desirable for cryptography, but it seems
unrelated to uses in gnuplot.

Consider the following script:

# Tabulate first 1000 numbers in sequences generated by successive seeds
set samples 1000
set table '1.dat'
print rand(101)
plot '+' using 0:(rand(0))
set table '2.dat'
print rand(102)
plot '+' using 0:(rand(0))
# Compare the paired elements within the two sequences
system("join 1.dat 2.dat > 1_2.dat")
plot "1_2.dat" using 2:4

No correlation is evident in the two sequences, even
though they are seeded with successive integers and thus
the numerical difference between their initial elements
is small. 

> 
> This is a problem if we want to use the standard trick of initializing
> the PRNG with the current time (as the user wanted to in [1]):
> gnuplot> print rand(real(system("date +%s")))
> 0.596925441003784
> gnuplot> print rand(real(system("date +%s")))
> 0.596924809567054
> gnuplot> print rand(real(system("date +%s")))
> 0.596924178130323
> gnuplot> print rand(real(system("date +%s")))
> 0.596923546693592

I don't see why this is a problem. Explain, please?
If the idea is to generate different streams of random numbers,
I think that it succeeds at that.
True, the initial values of sequences generated in close succession
are similar, but again I have to ask why this is a problem?
If it bothers you, can't you just start with the second value of each
sequence?   

> date +%s returns the time in the Unix time_t format (32 bit integer),
> however, the upper bits rarely change in that format.

So why not use system("date +%N")?


 
> 1) and 2) may be bugs in the implementation that are easy to fix, but 3)
> is a deeper problem. Of course, this kind of behavior can be expected
> from a linear congruence generator - but then maybe it's time to
> consider changing to a different algorithm.
> 
> Péter Juhász
> ps. Merry Christmas and a Happy New Year to you all!