Re: Help: Randomizing a List of Numbers

From: Tim Smith (reply_in_group_at_mouse-potato.com)
Date: 07/23/04 

Date: Thu, 22 Jul 2004 22:49:32 GMT

On 2004-07-22, Bill Emerson <no@one.home> wrote:
> How do I get the output of /dev/random into rn()? I really don't know how
> to write that function. Rand() won't do here.

Here's an implementation of rn(A,B) that takes two unsigned longs, A and B,
and returns an unsigned long in [A,B], for machines where unsigned long is
a 32 bit number. It uses /dev/urandom, but you can, of course, change it
to /dev/random. It reads 128*32 bits at a time from /dev/urandom, so as
to avoid the overhead of reading it each time you need a random number.
You'll have to stick in something for the two fatal error cases (can't open
/dev/urandom, and read fails).

#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <unistd.h>

unsigned long rn(unsigned long A, unsigned long B)
{
    static unsigned long bits[128];
    static int avail = 0;

    if ( avail == 0 )
    {
        int got;
        int fd = open( "/dev/urandom", O_RDONLY );
        if ( fd )
        {
            // FATAL ERROR!
        }
        got = read( fd, bits, sizeof bits );
        close( fd );
        if ( got != sizeof bits )
        {
            // FATAL ERROR!
        }
        avail = sizeof bits / sizeof bits[0];
    }

    return (unsigned long)(A + (B-A+1) * (bits[--avail] / 4294967296.0));
}

BTW, something to note. It would be easy to change this from taking
unsigned longs to taking longs or ints, and then you might think that you
could easily use it for negative numbers. E.g., rn(-10,10) to generate a
random integer in [-10,10]. And it would indeed do that...but the
distribution would be off. You'd actually get something in [-9,10] most of
the time. You'd get -10 only about 1/4294967296 of the time on most
computers.

This is because the people who design computers think that negative floats
should round toward zero when being converted to an integer, rather than the
more sensible rounding toward negative infinity, so to produce -10 in the
above example, bits[--avail] has to be exactly 0.

(This is also why one has to be careful of the % operator when negative
numbes are involved. You usually want something like -1%5 to be 4, but you
usually get -1. That comes from wanting to keep int(A/B)*B + A%B = A.
Since they round toward zero, when A<0, B>0, the int(A/B)*B ends up B higher
than it would be if they rounded toward negative infinity, and so the A%B
has to be in [-B+1,0] rather than [0,B-1] like you usually want). 

-- 
--Tim Smith