home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1992 #26 / NN_1992_26.iso / spool / comp / arch / 10397 < prev    next >
Encoding:
Text File  |  1992-11-04  |  13.2 KB  |  355 lines
  1. Path: sparky!uunet!ogicse!flop.ENGR.ORST.EDU!flop.ENGR.ORST.EDU!usenet
  2. From: crowl@jade.CS.ORST.EDU (Lawrence Crowl)
  3. Newsgroups: comp.arch
  4. Subject: Exponential Random Number Generator
  5. Message-ID: <1dam8vINNl5@flop.ENGR.ORST.EDU>
  6. Date: 5 Nov 92 08:33:35 GMT
  7. Article-I.D.: flop.1dam8vINNl5
  8. Distribution: world
  9. Organization: Computer Science Department, Oregon State University
  10. Lines: 342
  11. NNTP-Posting-Host: jade.cs.orst.edu
  12.  
  13. After spending substantial time with Herman Rubin's random number
  14. generator routines, which are based on acceptance-rejection methods, I
  15. reached the conclusion that they would never be as fast as more direct
  16. computation based on the inverse of the probability distribution
  17. functions.  I wrote the following program to generate random numbers
  18. using the latter method.  The routines are based primarily on inlined
  19. table-lookup of the inverse distribution.
  20.  
  21. I haven't tested the quality of the distributions, could someone who
  22. knows how please do so?
  23.  
  24. Compiled with gcc -O on a Sun SparcStation 2, this program produces
  25.  
  26.   1,317,523 uniform random numbers per second,
  27.     939,849 geometric random numbers per second,
  28.     975,609 poisson random numbers per second, or
  29.     199,401 exponential random numbers per second.
  30.  
  31. Some improvement can be gained with a find-first-bit-set operation, but
  32. it won't be a large gain.  I'm sure these routines have room for
  33. improvement, but they're pretty fast as it stands.
  34.  
  35. Note that the performance in these routines comes primarily from
  36. looking carefully at the problem and secondarily from inlining
  37. table-lookup.  Can we do better?
  38.  
  39.   Lawrence Crowl        503-737-2554    Computer Science Department
  40.                crowl@cs.orst.edu    Oregon State University
  41.           ...!hplabs!hp-pcd!orstcs!crowl    Corvallis, Oregon,  97331-3202
  42.  
  43. /*
  44.  
  45. Crowl_randoms.c
  46.  
  47. This program computes random numbers of various flavors
  48. based on the inverse probability distribution functions.
  49.  
  50. The intent is that these routines be fast,
  51. though there is still room for straightforward optimization.
  52.  
  53. This program assumes that unsigned ints are 32 bits.
  54.  
  55. In evalutating the work of the following routines,
  56. I assume that local variables will be in registers.
  57. With the exception of recursive calls, gcc's inlining is effective.
  58.  
  59. */
  60.  
  61. #include <stdio.h>
  62. #include <math.h>
  63.  
  64. /*
  65. Return a semi-uniform random unsigned integer variable
  66. modulo the representation of unsigned integers.
  67. Average work is one load, one store, one multiplication,
  68. one addition, and one return.
  69. */
  70.  
  71. static unsigned semiseed = 0 ;
  72.  
  73. inline unsigned semiuniform_unsigned( void )
  74.     {
  75.     return semiseed = (semiseed * 25173 + 13849) ;
  76.     }
  77.  
  78. /*
  79. Return a uniform random unsigned integer variable modulo the representation,
  80. where all bits are uniform and independent.
  81. I derived this routine from what Herman Rubin posted to comp.lang.misc.
  82. Average work for straightforward computation is nearly
  83. two comparisons, five loads, one store,
  84. three subtractions with one, one addition, and one return.
  85. */
  86.  
  87. #define SEED_WORDS 607
  88. #define SEED_TRAIL 460
  89. static unsigned seed_store[ SEED_WORDS ] ;
  90. static unsigned *current_rand ;
  91.     
  92. void reseed_uniform( void )
  93.     { /* this routine is vectorizable and inlineable */
  94.     unsigned *primary_seed = &(seed_store[ SEED_WORDS - 1 ]) ;
  95.     unsigned *secondary_seed = &(seed_store[ SEED_TRAIL - 1 ]) ;
  96.     while ( secondary_seed >= &(seed_store[ 0 ]) )
  97.         *primary_seed-- += *secondary_seed-- ;
  98.     secondary_seed = &(seed_store[ SEED_WORDS - 1 ]) ;
  99.     while ( primary_seed >= &(seed_store[ 0 ]) )
  100.         *primary_seed-- += *secondary_seed-- ;
  101.     current_rand = &(seed_store[ SEED_WORDS ]) ;
  102.     }
  103.  
  104. inline unsigned uniform_unsigned( void )
  105.     { /* this routine is easily inlined */
  106.     if ( current_rand != &(seed_store[ 0 ]) )
  107.         return *--current_rand ;
  108.     reseed_uniform( ) ;
  109.     return *(current_rand = &(seed_store[ SEED_WORDS - 1 ])) ;
  110.     }
  111.  
  112. void initialize_seed( void )
  113.     {
  114.     int i ;
  115.     for ( i = 0 ;  i < SEED_WORDS ;  i++ )
  116.         seed_store[ i ] = semiuniform_unsigned( ) ; /* use weak generator */
  117.     reseed_uniform( ) ; /* scramble after weak initialization */
  118.     }
  119.  
  120. /*
  121. Return a geometric random variable with mean 2.
  122. TeX: $p(y) = q^{y-1}p$ where $p=.5$ and $q=1-p$.
  123. It fails if the random variable exceeds the representation for unsigned,
  124. but then, you probably won't wait for it to fail.
  125. This routine relies on the independence of bits in uniform_unsigned.
  126. Average work is one call to uniform unsigned, one assignment,
  127. two comparisons, one return,
  128. and a smidgen of an add and a call to geometric_unsigned.
  129. Cliff Click suggested the recursion,
  130. rather than my earlier looping implementation.
  131. */
  132.  
  133. inline unsigned geometric_unsigned( void )
  134.     { /* this routine is easily inlined if a compiled version is available */
  135.     unsigned uniform = uniform_unsigned( ) ;
  136.     if ( uniform >= 0x80000000 ) return 1 ;
  137.     if ( uniform >= 0x40000000 ) return 2 ;
  138.     if ( uniform >= 0x20000000 ) return 3 ;
  139.     if ( uniform >= 0x10000000 ) return 4 ;
  140.     if ( uniform >= 0x08000000 ) return 5 ;
  141.     if ( uniform >= 0x04000000 ) return 6 ;
  142.     if ( uniform >= 0x02000000 ) return 7 ;
  143.     if ( uniform >= 0x01000000 ) return 8 ;
  144.     if ( uniform >= 0x00800000 ) return 9 ;
  145.     if ( uniform >= 0x00400000 ) return 10 ;
  146.     if ( uniform >= 0x00200000 ) return 11 ;
  147.     if ( uniform >= 0x00100000 ) return 12 ;
  148.     if ( uniform >= 0x00080000 ) return 13 ;
  149.     if ( uniform >= 0x00040000 ) return 14 ;
  150.     if ( uniform >= 0x00020000 ) return 15 ;
  151.     if ( uniform >= 0x00010000 ) return 16 ;
  152.     if ( uniform >= 0x00008000 ) return 17 ;
  153.     if ( uniform >= 0x00004000 ) return 18 ;
  154.     if ( uniform >= 0x00002000 ) return 19 ;
  155.     if ( uniform >= 0x00001000 ) return 20 ;
  156.     if ( uniform >= 0x00000800 ) return 21 ;
  157.     if ( uniform >= 0x00000400 ) return 22 ;
  158.     if ( uniform >= 0x00000200 ) return 23 ;
  159.     if ( uniform >= 0x00000100 ) return 24 ;
  160.     if ( uniform >= 0x00000080 ) return 25 ;
  161.     if ( uniform >= 0x00000040 ) return 26 ;
  162.     if ( uniform >= 0x00000020 ) return 27 ;
  163.     if ( uniform >= 0x00000010 ) return 28 ;
  164.     if ( uniform >= 0x00000008 ) return 29 ;
  165.     if ( uniform >= 0x00000004 ) return 30 ;
  166.     if ( uniform >= 0x00000002 ) return 31 ;
  167.     if ( uniform >= 0x00000001 ) return 32 ;
  168.     /* no set bit found (probability 1 in 2**32), go to next word */
  169.     return 32 + geometric_unsigned( ) ;
  170.     }
  171.  
  172. /*
  173. Return a Poisson random variable, with mean .5.
  174. TeX: $p(y) = {m^y\over y!}e^{-m}$ where $m=.5$.
  175. This routine is based on table lookup of
  176. the inverse of the probability distribution function.
  177. It interprets some unsigned numbers as fixed-point numbers in 0 <= x < 1.
  178. The average work is essentially:
  179. one call to uniform_unsigned,
  180. two comparisons with large constants,
  181. and one return.
  182. */
  183.  
  184. inline unsigned poisson_unsigned( void )
  185.     {
  186.     unsigned uniform = uniform_unsigned( ) ;
  187.     if ( uniform >= 0x9B4597E4 ) return  0 ;
  188.     if ( uniform >= 0x1368B2FD ) return  1 ;
  189.     if ( uniform >= 0x033C1DD5 ) return  2 ;
  190.     if ( uniform >= 0x006783BB ) return  3 ;
  191.     if ( uniform >= 0x000A59FA ) return  4 ;
  192.     if ( uniform >= 0x0000DCD5 ) return  5 ;
  193.     if ( uniform >= 0x00000FC7 ) return  6 ;
  194.     if ( uniform >= 0x000000FD ) return  7 ;
  195.     if ( uniform >= 0x0000000F ) return  8 ;
  196.     if ( uniform >= 0x00000001 ) return  9 ;
  197.     /* poisson not found with probability 1 in 2**32 */
  198.     /* get another uniform and keep trying */
  199.     /* think of the second uniform as another word of precision 
  200.        in a double-word uniform random variable */
  201.     uniform = uniform_unsigned( ) ;
  202.     if ( uniform >= 0xB3781489 ) return 10 ;
  203.     if ( uniform >= 0x08285E07 ) return 11 ;
  204.     if ( uniform >= 0x005703E9 ) return 12 ;
  205.     if ( uniform >= 0x000358C5 ) return 13 ;
  206.     if ( uniform >= 0x00001E9A ) return 14 ;
  207.     if ( uniform >= 0x00000106 ) return 15 ;
  208.     if ( uniform >= 0x00000009 ) return 16 ;
  209.     /* poisson not found with probability 9 in 2**64 */
  210.     /* get another uniform and keep trying */
  211.     uniform = uniform_unsigned( ) ;
  212.     if ( uniform >= 0x3D70045D ) return 17 ;
  213.     if ( uniform >= 0x01B4E3AE ) return 18 ;
  214.     if ( uniform >= 0x000B7F42 ) return 19 ;
  215.     if ( uniform >= 0x00004995 ) return 20 ;
  216.     if ( uniform >= 0x000001C1 ) return 21 ;
  217.     if ( uniform >= 0x0000000B ) return 22 ;
  218.     /* poisson not found with probability 11 in 2**96 */
  219.     /* get another uniform and keep trying */
  220.     uniform = uniform_unsigned( ) ;
  221.     if ( uniform >= 0x38BA0CF2 ) return 23 ;
  222.     if ( uniform >= 0x012E8AF0 ) return 24 ;
  223.     if ( uniform >= 0x00060D05 ) return 25 ;
  224.     if ( uniform >= 0x00001DCA ) return 26 ;
  225.     if ( uniform >= 0x0000008E ) return 27 ;
  226.     if ( uniform >= 0x00000003 ) return 28 ;
  227.     /* poisson not found with probability 3 in 2**128 */
  228.     /* we could get another uniform and keep going */
  229.     /* however, I quit and return the wrong answer */
  230.     return 29 ;
  231.     }
  232.  
  233. /*
  234. Return an exponential random double-precision floating point number,
  235. with a given mean.
  236. I derived this routine from Algorithm S on page 128 of
  237. Donald E. Knuth's "The Art of Computer Programming",
  238. volume 2 "Seminumerical Algorithms",
  239. second edition.
  240. It interprets some unsigned numbers as a fixed-point numbers in 0 <= x < 1.
  241.  
  242. This routine returns a double rather than the fixed point Herman Rubin wants.
  243. This is because fixed-point multiplication is a little hairy in C
  244. and the routine hasn't been tested.
  245.  
  246. I haven't evaluated the work for this routine,
  247. but I'd like to point out that on average
  248. it will compare fractional only twice.
  249. So, on average, this routine will be fast.
  250. */
  251.  
  252. inline double exponential_double( double mean )
  253.     {
  254.     unsigned newer, minimum ;
  255.     unsigned integral = geometric_unsigned( ) ;
  256.     unsigned fractional = uniform_unsigned( ) ;
  257.     if ( fractional < 0xB17217F8 /* (log(2)^1/1! */ )
  258.         return mean * (integral*M_LN2 + fractional/4294967296.0) ;
  259.     minimum = fractional ;
  260.     newer = uniform_unsigned( ) ;
  261.     if ( newer < minimum ) minimum = newer ;
  262.     if ( fractional < 0x3D7F7BFF /* (log(2)^2/2! */ )
  263.         return M_LN2 * mean * (integral + fractional/4294967296.0) ;
  264.     minimum = fractional ;
  265.     newer = uniform_unsigned( ) ;
  266.     if ( newer < minimum ) minimum = newer ;
  267.     if ( fractional < 0x0E35846C /* (log(2)^3/3! */ )
  268.         return M_LN2 * mean * (integral + fractional/4294967296.0) ;
  269.     minimum = fractional ;
  270.     newer = uniform_unsigned( ) ;
  271.     if ( newer < minimum ) minimum = newer ;
  272.     if ( fractional < 0x0276556E /* (log(2)^4/4! */ )
  273.         return M_LN2 * mean * (integral + fractional/4294967296.0) ;
  274.     minimum = fractional ;
  275.     newer = uniform_unsigned( ) ;
  276.     if ( newer < minimum ) minimum = newer ;
  277.     if ( fractional < 0x00576200 /* (log(2)^5/5! */ )
  278.         return M_LN2 * mean * (integral + fractional/4294967296.0) ;
  279.     minimum = fractional ;
  280.     newer = uniform_unsigned( ) ;
  281.     if ( newer < minimum ) minimum = newer ;
  282.     if ( fractional < 0x000A1849 /* (log(2)^6/6! */ )
  283.         return M_LN2 * mean * (integral + fractional/4294967296.0) ;
  284.     minimum = fractional ;
  285.     newer = uniform_unsigned( ) ;
  286.     if ( newer < minimum ) minimum = newer ;
  287.     if ( fractional < 0x0000FFE6 /* (log(2)^7/7! */ )
  288.         return M_LN2 * mean * (integral + fractional/4294967296.0) ;
  289.     minimum = fractional ;
  290.     newer = uniform_unsigned( ) ;
  291.     if ( newer < minimum ) minimum = newer ;
  292.     if ( fractional < 0x0000162C /* (log(2)^8/8! */ )
  293.         return M_LN2 * mean * (integral + fractional/4294967296.0) ;
  294.     minimum = fractional ;
  295.     newer = uniform_unsigned( ) ;
  296.     if ( newer < minimum ) minimum = newer ;
  297.     if ( fractional < 0x000001B5 /* (log(2)^9/9! */ )
  298.         return M_LN2 * mean * (integral + fractional/4294967296.0) ;
  299.     minimum = fractional ;
  300.     newer = uniform_unsigned( ) ;
  301.     if ( newer < minimum ) minimum = newer ;
  302.     if ( fractional < 0x0000001E /* (log(2)^10/10! */ )
  303.         return M_LN2 * mean * (integral + fractional/4294967296.0) ;
  304.     minimum = fractional ;
  305.     newer = uniform_unsigned( ) ;
  306.     if ( newer < minimum ) minimum = newer ;
  307.     if ( fractional < 0x00000002 /* (log(2)^11/11! */ )
  308.         return M_LN2 * mean * (integral + fractional/4294967296.0) ;
  309.     minimum = fractional ;
  310.     newer = uniform_unsigned( ) ;
  311.     if ( newer < minimum ) minimum = newer ;
  312.     /* terminate computation, error 2 parts in 2^32 */
  313.     return M_LN2 * mean * (integral + fractional/4294967296.0) ;
  314.     }
  315.  
  316. /*
  317. This is a non-testing main.
  318. */
  319.  
  320. int main( argc, argv )
  321.     int argc ;
  322.     char **argv ;
  323.     {
  324.     int count = 20 ;
  325.     char timing = ' ' ;
  326.     if ( argc > 1 ) count = atoi( argv[ 1 ] ) ;
  327.     if ( argc > 2 ) timing = argv[ 2 ][ 0 ] ;
  328.     initialize_seed( ) ;
  329.     if ( timing == 'u' )
  330.         while ( count-- > 0 ) uniform_unsigned( ) ;
  331.     else if ( timing == 'g' )
  332.         while ( count-- > 0 ) geometric_unsigned( ) ;
  333.     else if ( timing == 'p' )
  334.         while ( count-- > 0 ) poisson_unsigned( ) ;
  335.     else if ( timing == 'e' )
  336.         while ( count-- > 0 ) exponential_double( 2.0 ) ;
  337.     else
  338.         {
  339.         (void)printf( " uniform geometric poisson exponential\n" ) ;
  340.         while ( count-- > 0 )
  341.             {
  342.             (void)printf( "%08X %9u %7u %#11.8f\n",
  343.                           uniform_unsigned( ),
  344.                           geometric_unsigned( ),
  345.                           poisson_unsigned( ),
  346.                           exponential_double( 2.0 ) ) ;
  347.             }
  348.         }
  349.     return 0 ;
  350.     }
  351. -- 
  352.   Lawrence Crowl        503-737-2554    Computer Science Department
  353.                crowl@cs.orst.edu    Oregon State University
  354.           ...!hplabs!hp-pcd!orstcs!crowl    Corvallis, Oregon,  97331-3202
  355.