Random Numbers ...

[Lawrence Crowl]

Most of my need for random numbers comes in the form of a random integer modulo
some application specific value.  A lot of random number generators return
either a floating value, or some integer within the complete range of integers,
which I must then hack at to get it within the range of integers that I want.
Rather than do this all the time, I wrote a random number generator which
accepts the modulus and does the modulo calculations internally.  Here is the
result of that effort.  On a Sun 2/170, the 16-bit version takes 54 micro-
seconds and the 32-bit version takes 158 micro-seconds.


/*-------------------------------------------------------------------- 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.  Because twos complement
integers have more negative values than positive, this generator uses the
negative of the absolute value for some operations.

If the modulus is constant, an macro implementation of this generator will
eliminate a division.  In addition, the modulus parameter is evaluated only
once, so it will have no undue side effects.

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

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

/*----------------------------------------------------------------------- tools

Note that since the parameter to NEGABS is evaluated more than once, the
actual argument should have no side effects.
*/

#define NEGABS( i ) ((i) > 0 ? -(i) : (i))

/*-------------------------------------------------------------- 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 short seed16 ;

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

short rand16mod( modulus )
    short modulus ;
    {
    seed16 = (short)(seed16 * 25173) + 13849 ;
    return (short)NEGABS(seed16) / (short)(-32768 / modulus) ;
    }

/*-------------------------------------------------------------- 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 int seed32 ;

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

int rand32mod( modulus )
    int modulus ;
    {
    seed32 = (seed32 * 1103515245) + 13849 ;
    return NEGABS(seed32) / (-2147483648 / modulus) ;
    }

/*------------------------------------------------- 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 1000000

void main( )
    {
    int i, j ;
    struct timeval t1, t2, t3 ;
    long d1, d2 ;
    float di ;

    /* demonstrate 16-bit version */
    for ( i = 4 ;  i <= 8 ;  i++ )
        {
        for ( j = 0 ;  j <= 20 ;  j++ )
            (void)printf( " %d", rand16mod( i ) ) ;
        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( 32 ) ;
    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 = 4 ;  i <= 8 ;  i++ )
        {
        for ( j = 0 ;  j <= 20 ;  j++ )
            (void)printf( " %d", rand32mod( i ) ) ;
        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( 32 ) ;
    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