From: crowl@jade.CS.ORST.EDU (Lawrence Crowl)
Newsgroups: comp.arch
Subject: Exponential Random Number Generator
Date: 5 Nov 92 08:33:35 GMT
Article-I.D.: flop.1dam8vINNl5
Distribution: world
Organization: Computer Science Department, Oregon State University
NNTP-Posting-Host: jade.cs.orst.edu

After spending substantial time with Herman Rubin's random number
generator routines, which are based on acceptance-rejection methods, I
reached the conclusion that they would never be as fast as more direct
computation based on the inverse of the probability distribution
functions.  I wrote the following program to generate random numbers
using the latter method.  The routines are based primarily on inlined
table-lookup of the inverse distribution.

I haven't tested the quality of the distributions, could someone who
knows how please do so?

Compiled with gcc -O on a Sun SparcStation 2, this program produces

  1,317,523 uniform random numbers per second,
    939,849 geometric random numbers per second,
    975,609 poisson random numbers per second, or
    199,401 exponential random numbers per second.

Some improvement can be gained with a find-first-bit-set operation, but
it won't be a large gain.  I'm sure these routines have room for
improvement, but they're pretty fast as it stands.

Note that the performance in these routines comes primarily from
looking carefully at the problem and secondarily from inlining
table-lookup.  Can we do better?

  Lawrence Crowl		503-737-2554	Computer Science Department
			   crowl@cs.orst.edu	Oregon State University
	      ...!hplabs!hp-pcd!orstcs!crowl	Corvallis, Oregon,  97331-3202

/*

Crowl_randoms.c

This program computes random numbers of various flavors
based on the inverse probability distribution functions.

The intent is that these routines be fast,
though there is still room for straightforward optimization.

This program assumes that unsigned ints are 32 bits.

In evalutating the work of the following routines,
I assume that local variables will be in registers.
With the exception of recursive calls, gcc's inlining is effective.

*/

#include <stdio.h>
#include <math.h>

/*
Return a semi-uniform random unsigned integer variable
modulo the representation of unsigned integers.
Average work is one load, one store, one multiplication,
one addition, and one return.
*/

static unsigned semiseed = 0 ;

inline unsigned semiuniform_unsigned( void )
    {
    return semiseed = (semiseed * 25173 + 13849) ;
    }

/*
Return a uniform random unsigned integer variable modulo the representation,
where all bits are uniform and independent.
I derived this routine from what Herman Rubin posted to comp.lang.misc.
Average work for straightforward computation is nearly
two comparisons, five loads, one store,
three subtractions with one, one addition, and one return.
*/

#define SEED_WORDS 607
#define SEED_TRAIL 460
static unsigned seed_store[ SEED_WORDS ] ;
static unsigned *current_rand ;
    
void reseed_uniform( void )
    { /* this routine is vectorizable and inlineable */
    unsigned *primary_seed = &(seed_store[ SEED_WORDS - 1 ]) ;
    unsigned *secondary_seed = &(seed_store[ SEED_TRAIL - 1 ]) ;
    while ( secondary_seed >= &(seed_store[ 0 ]) )
        *primary_seed-- += *secondary_seed-- ;
    secondary_seed = &(seed_store[ SEED_WORDS - 1 ]) ;
    while ( primary_seed >= &(seed_store[ 0 ]) )
        *primary_seed-- += *secondary_seed-- ;
    current_rand = &(seed_store[ SEED_WORDS ]) ;
    }

inline unsigned uniform_unsigned( void )
    { /* this routine is easily inlined */
    if ( current_rand != &(seed_store[ 0 ]) )
        return *--current_rand ;
    reseed_uniform( ) ;
    return *(current_rand = &(seed_store[ SEED_WORDS - 1 ])) ;
    }

void initialize_seed( void )
    {
    int i ;
    for ( i = 0 ;  i < SEED_WORDS ;  i++ )
        seed_store[ i ] = semiuniform_unsigned( ) ; /* use weak generator */
    reseed_uniform( ) ; /* scramble after weak initialization */
    }

/*
Return a geometric random variable with mean 2.
TeX: $p(y) = q^{y-1}p$ where $p=.5$ and $q=1-p$.
It fails if the random variable exceeds the representation for unsigned,
but then, you probably won't wait for it to fail.
This routine relies on the independence of bits in uniform_unsigned.
Average work is one call to uniform unsigned, one assignment,
two comparisons, one return,
and a smidgen of an add and a call to geometric_unsigned.
Cliff Click suggested the recursion,
rather than my earlier looping implementation.
*/

inline unsigned geometric_unsigned( void )
    { /* this routine is easily inlined if a compiled version is available */
    unsigned uniform = uniform_unsigned( ) ;
    if ( uniform >= 0x80000000 ) return 1 ;
    if ( uniform >= 0x40000000 ) return 2 ;
    if ( uniform >= 0x20000000 ) return 3 ;
    if ( uniform >= 0x10000000 ) return 4 ;
    if ( uniform >= 0x08000000 ) return 5 ;
    if ( uniform >= 0x04000000 ) return 6 ;
    if ( uniform >= 0x02000000 ) return 7 ;
    if ( uniform >= 0x01000000 ) return 8 ;
    if ( uniform >= 0x00800000 ) return 9 ;
    if ( uniform >= 0x00400000 ) return 10 ;
    if ( uniform >= 0x00200000 ) return 11 ;
    if ( uniform >= 0x00100000 ) return 12 ;
    if ( uniform >= 0x00080000 ) return 13 ;
    if ( uniform >= 0x00040000 ) return 14 ;
    if ( uniform >= 0x00020000 ) return 15 ;
    if ( uniform >= 0x00010000 ) return 16 ;
    if ( uniform >= 0x00008000 ) return 17 ;
    if ( uniform >= 0x00004000 ) return 18 ;
    if ( uniform >= 0x00002000 ) return 19 ;
    if ( uniform >= 0x00001000 ) return 20 ;
    if ( uniform >= 0x00000800 ) return 21 ;
    if ( uniform >= 0x00000400 ) return 22 ;
    if ( uniform >= 0x00000200 ) return 23 ;
    if ( uniform >= 0x00000100 ) return 24 ;
    if ( uniform >= 0x00000080 ) return 25 ;
    if ( uniform >= 0x00000040 ) return 26 ;
    if ( uniform >= 0x00000020 ) return 27 ;
    if ( uniform >= 0x00000010 ) return 28 ;
    if ( uniform >= 0x00000008 ) return 29 ;
    if ( uniform >= 0x00000004 ) return 30 ;
    if ( uniform >= 0x00000002 ) return 31 ;
    if ( uniform >= 0x00000001 ) return 32 ;
    /* no set bit found (probability 1 in 2**32), go to next word */
    return 32 + geometric_unsigned( ) ;
    }

/*
Return a Poisson random variable, with mean .5.
TeX: $p(y) = {m^y\over y!}e^{-m}$ where $m=.5$.
This routine is based on table lookup of
the inverse of the probability distribution function.
It interprets some unsigned numbers as fixed-point numbers in 0 <= x < 1.
The average work is essentially:
one call to uniform_unsigned,
two comparisons with large constants,
and one return.
*/

inline unsigned poisson_unsigned( void )
    {
    unsigned uniform = uniform_unsigned( ) ;
    if ( uniform >= 0x9B4597E4 ) return  0 ;
    if ( uniform >= 0x1368B2FD ) return  1 ;
    if ( uniform >= 0x033C1DD5 ) return  2 ;
    if ( uniform >= 0x006783BB ) return  3 ;
    if ( uniform >= 0x000A59FA ) return  4 ;
    if ( uniform >= 0x0000DCD5 ) return  5 ;
    if ( uniform >= 0x00000FC7 ) return  6 ;
    if ( uniform >= 0x000000FD ) return  7 ;
    if ( uniform >= 0x0000000F ) return  8 ;
    if ( uniform >= 0x00000001 ) return  9 ;
    /* poisson not found with probability 1 in 2**32 */
    /* get another uniform and keep trying */
    /* think of the second uniform as another word of precision 
       in a double-word uniform random variable */
    uniform = uniform_unsigned( ) ;
    if ( uniform >= 0xB3781489 ) return 10 ;
    if ( uniform >= 0x08285E07 ) return 11 ;
    if ( uniform >= 0x005703E9 ) return 12 ;
    if ( uniform >= 0x000358C5 ) return 13 ;
    if ( uniform >= 0x00001E9A ) return 14 ;
    if ( uniform >= 0x00000106 ) return 15 ;
    if ( uniform >= 0x00000009 ) return 16 ;
    /* poisson not found with probability 9 in 2**64 */
    /* get another uniform and keep trying */
    uniform = uniform_unsigned( ) ;
    if ( uniform >= 0x3D70045D ) return 17 ;
    if ( uniform >= 0x01B4E3AE ) return 18 ;
    if ( uniform >= 0x000B7F42 ) return 19 ;
    if ( uniform >= 0x00004995 ) return 20 ;
    if ( uniform >= 0x000001C1 ) return 21 ;
    if ( uniform >= 0x0000000B ) return 22 ;
    /* poisson not found with probability 11 in 2**96 */
    /* get another uniform and keep trying */
    uniform = uniform_unsigned( ) ;
    if ( uniform >= 0x38BA0CF2 ) return 23 ;
    if ( uniform >= 0x012E8AF0 ) return 24 ;
    if ( uniform >= 0x00060D05 ) return 25 ;
    if ( uniform >= 0x00001DCA ) return 26 ;
    if ( uniform >= 0x0000008E ) return 27 ;
    if ( uniform >= 0x00000003 ) return 28 ;
    /* poisson not found with probability 3 in 2**128 */
    /* we could get another uniform and keep going */
    /* however, I quit and return the wrong answer */
    return 29 ;
    }

/*
Return an exponential random double-precision floating point number,
with a given mean.
I derived this routine from Algorithm S on page 128 of
Donald E. Knuth's "The Art of Computer Programming",
volume 2 "Seminumerical Algorithms",
second edition.
It interprets some unsigned numbers as a fixed-point numbers in 0 <= x < 1.

This routine returns a double rather than the fixed point Herman Rubin wants.
This is because fixed-point multiplication is a little hairy in C
and the routine hasn't been tested.

I haven't evaluated the work for this routine,
but I'd like to point out that on average
it will compare fractional only twice.
So, on average, this routine will be fast.
*/

inline double exponential_double( double mean )
    {
    unsigned newer, minimum ;
    unsigned integral = geometric_unsigned( ) ;
    unsigned fractional = uniform_unsigned( ) ;
    if ( fractional < 0xB17217F8 /* (log(2)^1/1! */ )
        return mean * (integral*M_LN2 + fractional/4294967296.0) ;
    minimum = fractional ;
    newer = uniform_unsigned( ) ;
    if ( newer < minimum ) minimum = newer ;
    if ( fractional < 0x3D7F7BFF /* (log(2)^2/2! */ )
        return M_LN2 * mean * (integral + fractional/4294967296.0) ;
    minimum = fractional ;
    newer = uniform_unsigned( ) ;
    if ( newer < minimum ) minimum = newer ;
    if ( fractional < 0x0E35846C /* (log(2)^3/3! */ )
        return M_LN2 * mean * (integral + fractional/4294967296.0) ;
    minimum = fractional ;
    newer = uniform_unsigned( ) ;
    if ( newer < minimum ) minimum = newer ;
    if ( fractional < 0x0276556E /* (log(2)^4/4! */ )
        return M_LN2 * mean * (integral + fractional/4294967296.0) ;
    minimum = fractional ;
    newer = uniform_unsigned( ) ;
    if ( newer < minimum ) minimum = newer ;
    if ( fractional < 0x00576200 /* (log(2)^5/5! */ )
        return M_LN2 * mean * (integral + fractional/4294967296.0) ;
    minimum = fractional ;
    newer = uniform_unsigned( ) ;
    if ( newer < minimum ) minimum = newer ;
    if ( fractional < 0x000A1849 /* (log(2)^6/6! */ )
        return M_LN2 * mean * (integral + fractional/4294967296.0) ;
    minimum = fractional ;
    newer = uniform_unsigned( ) ;
    if ( newer < minimum ) minimum = newer ;
    if ( fractional < 0x0000FFE6 /* (log(2)^7/7! */ )
        return M_LN2 * mean * (integral + fractional/4294967296.0) ;
    minimum = fractional ;
    newer = uniform_unsigned( ) ;
    if ( newer < minimum ) minimum = newer ;
    if ( fractional < 0x0000162C /* (log(2)^8/8! */ )
        return M_LN2 * mean * (integral + fractional/4294967296.0) ;
    minimum = fractional ;
    newer = uniform_unsigned( ) ;
    if ( newer < minimum ) minimum = newer ;
    if ( fractional < 0x000001B5 /* (log(2)^9/9! */ )
        return M_LN2 * mean * (integral + fractional/4294967296.0) ;
    minimum = fractional ;
    newer = uniform_unsigned( ) ;
    if ( newer < minimum ) minimum = newer ;
    if ( fractional < 0x0000001E /* (log(2)^10/10! */ )
        return M_LN2 * mean * (integral + fractional/4294967296.0) ;
    minimum = fractional ;
    newer = uniform_unsigned( ) ;
    if ( newer < minimum ) minimum = newer ;
    if ( fractional < 0x00000002 /* (log(2)^11/11! */ )
        return M_LN2 * mean * (integral + fractional/4294967296.0) ;
    minimum = fractional ;
    newer = uniform_unsigned( ) ;
    if ( newer < minimum ) minimum = newer ;
    /* terminate computation, error 2 parts in 2^32 */
    return M_LN2 * mean * (integral + fractional/4294967296.0) ;
    }

/*
This is a non-testing main.
*/

int main( argc, argv )
    int argc ;
    char **argv ;
    {
    int count = 20 ;
    char timing = ' ' ;
    if ( argc > 1 ) count = atoi( argv[ 1 ] ) ;
    if ( argc > 2 ) timing = argv[ 2 ][ 0 ] ;
    initialize_seed( ) ;
    if ( timing == 'u' )
        while ( count-- > 0 ) uniform_unsigned( ) ;
    else if ( timing == 'g' )
        while ( count-- > 0 ) geometric_unsigned( ) ;
    else if ( timing == 'p' )
        while ( count-- > 0 ) poisson_unsigned( ) ;
    else if ( timing == 'e' )
        while ( count-- > 0 ) exponential_double( 2.0 ) ;
    else
        {
        (void)printf( " uniform geometric poisson exponential\n" ) ;
        while ( count-- > 0 )
            {
            (void)printf( "%08X %9u %7u %#11.8f\n",
                          uniform_unsigned( ),
                          geometric_unsigned( ),
                          poisson_unsigned( ),
                          exponential_double( 2.0 ) ) ;
            }
        }
    return 0 ;
    }
-- 
  Lawrence Crowl		503-737-2554	Computer Science Department
			   crowl@cs.orst.edu	Oregon State University
	      ...!hplabs!hp-pcd!orstcs!crowl	Corvallis, Oregon,  97331-3202
Newsgroups: comp.arch