Path: utzoo!mnetor!uunet!lll-winken!lll-lcc!ames!ll-xn!husc6!bbn!rochester!crowl From: crowl@cs.rochester.edu (Lawrence Crowl) Newsgroups: comp.lang.c Subject: Re: Random Numbers ... Message-ID: <7135@sol.ARPA> Date: 25 Feb 88 22:39:33 GMT References: <11972@brl-adm.ARPA> <7097@sol.ARPA> <5555@cit-vax.Caltech.Edu> Reply-To: crowl@cs.rochester.edu (Lawrence Crowl) Organization: U of Rochester, CS Dept, Rochester, NY Lines: 188 In article <5555@cit-vax.Caltech.Edu> wen-king@cit-vlsi.UUCP (Wen-King Su) writes: >In article <7097@sol.ARPA> crowl@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@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 #include #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@cs.rochester.edu Computer Science Department ...!{allegra,decvax,rutgers}!rochester!crowl Rochester, New York, 14627