Random Numbers ...

[Lawrence Crowl]

In article <5555 at cit-vax.Caltech.Edu> wen-king at cit-vlsi.UUCP
(Wen-King Su) writes:
>In article <7097 at sol.ARPA> crowl at cs.rochester.edu (Lawrence Crowl) writes:
><#define NEGABS( i ) ((i) > 0 ? -(i) : (i))
>>short rand16mod( modulus )
><    short modulus ;
>>    {
><    seed16 = (short)(seed16 * 25173) + 13849 ;
>>    return (short)NEGABS(seed16) / (short)(-32768 / modulus) ;
>            ^^^^^^^^^^^^^^^^^^^^^
>Don't do this.  Now you have a random number source that does not have
>an uniform distribution; both '0' and MAX_NEGATIVE occur at half the
>frequency as the other numbers.  Use unsigned integers instead.

Since this was developed for a language providing only signed integers, the
use of unsigned numbers did not occur to me.  The asymmetry with reguards to
0 and MAX_NEGATIVE are lost in the noise since the codes purpose was speed,
not extreme quality.  The use of MAX_NEGATIVE provides the largest signed
power of two.  I have fixed this since it costs little.  (Given the loss of
inlining brought about by the loop introduced to fix the following problem.)

>Even if the random source is uniform over a range, say R, you can't expect
>the the result of a simple division to yield another uniform distribution
>unless R is a multiple of the divisor.  For example, suppose the random
>source has a range of 0-7, and the divisor happens to be 3 (so that you get
>the numbers 0, 1, and 2 randomly).

If you look carefully at my old code, you will find that it uniformly
generates numbers in the range 0 .. modulus-1.  Unfortunately, it also
occasionally generates the value modulus.  This is not just a distribution
problem, it is an outright bug.

Taking Wen-King Su suggestions, I have redone the previous posting and
provide the result here.  Please throw out the old version and forget I ever
posted it!  The new code takes 60 micro-seconds for the 16-bit version and
150 micro-seconds for the 32-bit version on a Sun 2/170.


/*-------------------------------------------------------------------- random.c

This program contains a linear congruential pseudo-random number generator.
Each random number generation takes a parameter, modulus, and returns a
pseudo-random number in the range 0 ... modulus-1.  The pseudo-random number
is taken from the high-order bits of the seed.

This file contains two versions of the generator, one for 16-bit arithmetic,
and one for 32-bit arithmetic.

Lawrence Crowl
University of Rochester
Computer Science Department
Rochester, New York, 14627
*/

/*-------------------------------------------------------------- 16-bit version

This version assumes that shorts are 16 bits.  This allows the code to run
on machines with 32-bit ints but 16-bit hardware, such as the 68000.  The
extensive casting allows reasonable compilers to generate 16-bit multiply
instructions.

The coefficients for the generator are from Peter Grogono, "Programming in
Pascal", Section 4.5, Addison-Wesley Publishing Company, 1978.  I do not
know if these coefficients have been tested.
*/

static unsigned short seed16 ;

void rand16seed( seed )
    unsigned short seed ;
    {
    seed16 = seed ;
    }

short rand16mod( modulus )
    register unsigned short modulus ;
    {
    register unsigned short result, reg_seed ;

    reg_seed = seed16 ;
    do  {
        reg_seed = (unsigned short)(reg_seed * 25173) + 13849 ;
        result = (unsigned short)(reg_seed >> 1)
               / (unsigned short)(32767 / modulus) ;
        /* result = reg_seed / (unsigned short)(65535 / modulus) ; */
        /* the Sun compiler thinks 65535 needs long division */
        }
    while ( result >= modulus ) ;
    seed16 = reg_seed ;
    return result ;
    }

/*-------------------------------------------------------------- 32-bit version

This version assumes that ints are 32 bits.

The coefficients for the generator are from David Dantowitz in article
<11972 at brl-adm.ARPA> of the comp.lang.c news group.  I do not know if
these coefficients have been tested.
*/

static unsigned int seed32 ;

void rand32seed( seed )
    unsigned int seed ;
    {
    seed32 = seed ;
    }

int rand32mod( modulus )
    register unsigned int modulus ;
    {
    register unsigned int result, reg_seed ;

    reg_seed = seed32 ;
    do  {
        reg_seed = (reg_seed * 1103515245) + 12345 ;
        result = reg_seed / (4294967295 / modulus) ;
        }
    while ( result >= modulus ) ;
    seed32 = reg_seed ;
    return result ;
    }

/*------------------------------------------------- demonstration program

This is not a test program.  For test procedures, see Donald E. Knuth,
"The Art of Computer Programming", Chaper 3, Volume 2, Second Edition,
Addison-Wesley Publishing Company, 1981
*/

#include <stdio.h>
#include <sys/time.h>

#define ITERATIONS 100000

void main( argc, argv )
    int argc ;
    char **argv ;
    {
    int i ;
    struct timeval t1, t2, t3 ;
    long d1, d2 ;
    float di ;
    int samples = atoi( argv[ 1 ] ) ;
    int modulus = atoi( argv[ 2 ] ) ;

    /* demonstrate 16-bit version */
    for ( i = 0 ;  i < samples ;  i++ )
        (void)printf( " %d", rand16mod( modulus ) ) ;
    printf( "\n" ) ;

    /* time 16-bit version */
    gettimeofday( &t1, (struct timezone *)NULL ) ;
    for ( i = 0 ;  i < ITERATIONS ;  i++ )
        ;
    gettimeofday( &t2, (struct timezone *)NULL ) ;
    for ( i = 0 ;  i < ITERATIONS ;  i++ )
        (void)rand16mod( modulus ) ;
    gettimeofday( &t3, (struct timezone *)NULL ) ;
    d1 = (t2.tv_sec - t1.tv_sec) * 1000000 + (t2.tv_usec - t1.tv_usec) ;
    d2 = (t3.tv_sec - t2.tv_sec) * 1000000 + (t3.tv_usec - t2.tv_usec) ;
    di = (d2 - d1) / (float) ITERATIONS ;
    (void)printf( "microseconds 16-bit, null %d body %d iteration %f\n",
                  d1, d2, di ) ;

    /* demonstrate 32-bit version */
    for ( i = 0 ;  i < samples ;  i++ )
        (void)printf( " %d", rand32mod( modulus ) ) ;
    printf( "\n" ) ;

    /* time 32-bit version */
    gettimeofday( &t1, (struct timezone *)NULL ) ;
    for ( i = 0 ;  i < ITERATIONS ;  i++ )
        ;
    gettimeofday( &t2, (struct timezone *)NULL ) ;
    for ( i = 0 ;  i < ITERATIONS ;  i++ )
        (void)rand32mod( modulus ) ;
    gettimeofday( &t3, (struct timezone *)NULL ) ;
    d1 = (t2.tv_sec - t1.tv_sec) * 1000000 + (t2.tv_usec - t1.tv_usec) ;
    d2 = (t3.tv_sec - t2.tv_sec) * 1000000 + (t3.tv_usec - t2.tv_usec) ;
    di = (d2 - d1) / (float) ITERATIONS ;
    (void)printf( "microseconds 32-bit, null %d body %d iteration %f\n",
                  d1, d2, di ) ;
    }
-- 
  Lawrence Crowl		716-275-9499	University of Rochester
		      crowl at cs.rochester.edu	Computer Science Department
...!{allegra,decvax,rutgers}!rochester!crowl	Rochester, New York,  14627