home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS - Coast to Coast / simteldosarchivecoasttocoast2.iso / biology / gsrc208a.zip / R250.C < prev    next >
C/C++ Source or Header  |  1992-12-05  |  7KB  |  195 lines
  1. /******************************************************************************
  2. *  Module:  r250.c
  3. *  Description: implements R250 random number generator, from S.
  4. *  Kirkpatrick and E. Stoll, Journal of Computational Physics, 40,
  5. *  p. 517 (1981).
  6. *  Written by:    W. L. Maier
  7. *
  8. *    METHOD...
  9. *        16 parallel copies of a linear shift register with
  10. *        period 2^250 - 1.  FAR longer period than the usual
  11. *        linear congruent generator, and commonly faster as
  12. *        well.  (For details see the above paper, and the
  13. *        article in DDJ referenced below.)
  14. *
  15. *    HISTORY...
  16. *        Sep 92: Number returned by dr250() is in range [0,1) instead
  17. *            of [0,1], so for example a random angle in the
  18. *            interval [0, 2*PI) can be calculated
  19. *            conveniently.  (J. R. Van Zandt <jrv@mbunix.mitre.org>)
  20. *        Aug 92: Initialization is optional.  Default condition is
  21. *            equivalent to initializing with the seed 12345,
  22. *            so that the sequence of random numbers begins:
  23. *            1173, 53403, 52352, 35341...  (9996 numbers
  24. *            skipped) ...57769, 14511, 46930, 11942, 7978,
  25. *            56163, 46506, 45768, 21162, 43113...  Using ^=
  26. *            operator rather than two separate statements.
  27. *            Initializing with own linear congruent
  28. *            pseudorandom number generator for portability.
  29. *            Function prototypes moved to a header file.
  30. *            Implemented r250n() to generate numbers
  31. *            uniformly distributed in a specific range
  32. *            [0,n), because the common expedient rand()%n is
  33. *            incorrect.  (J. R. Van Zandt <jrv@mbunix.mitre.org>)
  34. *        May 91: Published by W. L. Maier, "A Fast Pseudo Random Number
  35. *            Generator", Dr. Dobb's Journal #176.
  36. ******************************************************************************/
  37.  
  38. #include <stdlib.h>
  39. #include "r250.h"
  40.  
  41. /**** Static variables ****/
  42. static int r250_index = 0;
  43. static unsigned int r250_buffer[250] = {
  44.     15301,57764,10921,56345,19316,43154,54727,49252,32360,49582,
  45.     26124,25833,34404,11030,26232,13965,16051,63635,55860,5184,
  46.     15931,39782,16845,11371,38624,10328,9139,1684,48668,59388,
  47.     13297,1364,56028,15687,63279,27771,5277,44628,31973,46977,
  48.     16327,23408,36065,52272,33610,61549,58364,3472,21367,56357,
  49.     56345,54035,7712,55884,39774,10241,50164,47995,1718,46887,
  50.     47892,6010,29575,54972,30458,21966,54449,10387,4492,644,
  51.     57031,41607,61820,54588,40849,54052,59875,43128,50370,44691,
  52.     286,12071,3574,61384,15592,45677,9711,23022,35256,45493,
  53.     48913,146,9053,5881,36635,43280,53464,8529,34344,64955,
  54.     38266,12730,101,16208,12607,58921,22036,8221,31337,11984,
  55.     20290,26734,19552,48,31940,43448,34762,53344,60664,12809,
  56.     57318,17436,44730,19375,30,17425,14117,5416,23853,55783,
  57.     57995,32074,26526,2192,11447,11,53446,35152,64610,64883,
  58.     26899,25357,7667,3577,39414,51161,4,58427,57342,58557,
  59.     53233,1066,29237,36808,19370,17493,37568,3,61468,38876,
  60.     17586,64937,21716,56472,58160,44955,55221,63880,1,32200,
  61.     62066,22911,24090,10438,40783,36364,14999,2489,43284,9898,
  62.     39612,9245,593,34857,41054,30162,65497,53340,27209,45417,
  63.     37497,4612,58397,52910,56313,62716,22377,40310,15190,34471,
  64.     64005,18090,11326,50839,62901,59284,5580,15231,9467,13161,
  65.     58500,7259,317,50968,2962,23006,32280,6994,18751,5148,
  66.     52739,49370,51892,18552,52264,54031,2804,17360,1919,19639,
  67.     2323,9448,43821,11022,45500,31509,49180,35598,38883,19754,
  68.     987,11521,55494,38056,20664,2629,50986,31009,54043,59743
  69.     };
  70.  
  71. static unsigned myrand();
  72. static void mysrand(unsigned newseed);
  73.  
  74. /**** Function: r250_init
  75.     Description: initializes r250 random number generator. ****/
  76. void r250_init(int seed)
  77. {
  78. /*---------------------------------------------------------------------------*/
  79.     int        j, k;
  80.     unsigned int mask;
  81.     unsigned int msb;
  82. /*---------------------------------------------------------------------------*/
  83.     mysrand(seed);
  84.     r250_index = 0;
  85.     for (j = 0; j < 250; j++)     /* Fill the r250 buffer with 15-bit values */
  86.         r250_buffer[j] = myrand();
  87.     for (j = 0; j < 250; j++)     /* Set some of the MS bits to 1 */
  88.         if (myrand() > 16384)
  89.             r250_buffer[j] |= 0x8000;
  90.     msb = 0x8000;       /* To turn on the diagonal bit   */
  91.     mask = 0xffff;      /* To turn off the leftmost bits */
  92.     for (j = 0; j < 16; j++)
  93.         {
  94.         k = 11 * j + 3;             /* Select a word to operate on        */
  95.         r250_buffer[k] &= mask;     /* Turn off bits left of the diagonal */
  96.         r250_buffer[k] |= msb;      /* Turn on the diagonal bit           */
  97.         mask >>= 1;
  98.         msb >>= 1;
  99.         }
  100. }
  101.  
  102. /**** Function: r250 Description: returns a random unsigned integer k
  103.                 uniformly distributed in the interval 0 <= k < 65536.  ****/
  104. unsigned int r250()
  105. {
  106. /*---------------------------------------------------------------------------*/
  107.     register int    j;
  108.     register unsigned int new_rand;
  109. /*---------------------------------------------------------------------------*/
  110.     if (r250_index >= 147)
  111.         j = r250_index - 147;      /* Wrap pointer around */
  112.     else
  113.         j = r250_index + 103;
  114.  
  115.     new_rand = r250_buffer[r250_index] ^= r250_buffer[j];
  116.  
  117.     if (r250_index >= 249)      /* Increment pointer for next time */
  118.         r250_index = 0;
  119.     else
  120.         r250_index++;
  121.  
  122.     return new_rand;
  123. }
  124.  
  125. /**** Function: r250n Description: returns a random unsigned integer k
  126.                     uniformly distributed in the interval 0 <= k < n ****/
  127. unsigned int r250n(unsigned n)
  128. {
  129. /*---------------------------------------------------------------------------*/
  130.     register int    j;
  131.     register unsigned int new_rand, limit;
  132. /*---------------------------------------------------------------------------*/
  133.     limit = (65535U/n)*n;
  134.     do
  135.         {new_rand = r250();
  136.         if (r250_index >= 147)
  137.             j = r250_index - 147;      /* Wrap pointer around */
  138.         else
  139.             j = r250_index + 103;
  140.  
  141.         new_rand = r250_buffer[r250_index] ^= r250_buffer[j];
  142.  
  143.         if (r250_index >= 249)      /* Increment pointer for next time */
  144.             r250_index = 0;
  145.         else
  146.             r250_index++;
  147.         } while(new_rand >= limit);
  148.     return new_rand%n;
  149. }
  150.  
  151. /**** Function: dr250
  152.         Description: returns a random double z in range 0 <= z < 1.  ****/
  153. double dr250()
  154. {
  155. /*---------------------------------------------------------------------------*/
  156.     register int    j;
  157.     register unsigned int new_rand;
  158. /*---------------------------------------------------------------------------*/
  159.     if (r250_index >= 147)
  160.         j = r250_index - 147;     /* Wrap pointer around */
  161.     else
  162.         j = r250_index + 103;
  163.  
  164.     new_rand = r250_buffer[r250_index] ^= r250_buffer[j];
  165.  
  166.     if (r250_index >= 249)      /* Increment pointer for next time */
  167.         r250_index = 0;
  168.     else
  169.         r250_index++;
  170.     return new_rand / 65536.;   /* Return a number in [0.0 to 1.0) */
  171. }
  172.  
  173.  
  174. /***  linear congruent pseudorandom number generator for initialization ***/
  175.  
  176. static unsigned long seed=1;
  177.  
  178.     /*    Return a pseudorandom number in the interval 0 <= n < 32768.
  179.         This produces the following sequence of pseudorandom
  180.         numbers:
  181.         346, 130, 10982, 1090...  (9996 numbers skipped) ...23369,
  182.         2020, 5703, 12762, 10828, 16252, 28648, 27041, 23444, 6604...
  183.     */
  184. static unsigned myrand()
  185. {
  186.     seed = seed*0x015a4e35L + 1;
  187.     return (seed>>16)&0x7fff;
  188. }
  189.  
  190.     /*    Initialize above linear congruent pseudorandom number generator */
  191. static void mysrand(unsigned newseed)
  192. {    seed = newseed;
  193. }
  194.  
  195.