home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS - Coast to Coast / simteldosarchivecoasttocoast2.iso / ddjmag / ddj9105.zip / PSEUDO.ASC < prev    next >
Text File  |  1991-03-15  |  5KB  |  157 lines
  1. _A FAST PSEUDO RANDOM NUMBER GENERATOR_
  2. by W.L. Maier
  3.  
  4. [LISTING ONE]
  5.  
  6. /******************************************************************************
  7. *  Module:  r250.cpp   Description: implements R250 random number generator, 
  8. *  from S. Kirkpatrick and E. Stoll, Journal of Computational Physics, 40,
  9. *  p. 517 (1981). Written by:    W. L. Maier
  10. ******************************************************************************/
  11.  
  12. #include <stdlib.h>
  13.  
  14. /**** Static variables ****/
  15. static unsigned int r250_buffer[250];
  16. static int r250_index;
  17.  
  18. /**** Function prototypes  ****/
  19. void r250_init(int seed);
  20. unsigned int r250();
  21. double dr250();
  22.  
  23. /**** Function: r250_init  Description: initializes r250 random number 
  24. generator. ****/
  25. void r250_init(int seed)
  26. {
  27. /*---------------------------------------------------------------------------*/
  28.     int        j, k;
  29.     unsigned int mask;
  30.     unsigned int msb;
  31. /*---------------------------------------------------------------------------*/
  32.     srand(seed);
  33.     r250_index = 0;
  34.     for (j = 0; j < 250; j++)     /* Fill the r250 buffer with 15-bit values */
  35.         r250_buffer[j] = rand();
  36.     for (j = 0; j < 250; j++)     /* Set some of the MS bits to 1 */
  37.         if (rand() > 16384)
  38.             r250_buffer[j] |= 0x8000;
  39.     msb = 0x8000;       /* To turn on the diagonal bit   */
  40.     mask = 0xffff;      /* To turn off the leftmost bits */
  41.     for (j = 0; j < 16; j++)
  42.         {
  43.         k = 11 * j + 3;             /* Select a word to operate on        */
  44.         r250_buffer[k] &= mask;     /* Turn off bits left of the diagonal */
  45.         r250_buffer[k] |= msb;      /* Turn on the diagonal bit           */
  46.         mask >>= 1;
  47.         msb >>= 1;
  48.         }
  49. }
  50. /**** Function:  r250  Description: returns a random unsigned integer. ****/
  51. unsigned int r250()
  52. {
  53. /*---------------------------------------------------------------------------*/
  54.     register int    j;
  55.     register unsigned int new_rand;
  56. /*---------------------------------------------------------------------------*/
  57.     if (r250_index >= 147)
  58.         j = r250_index - 147;      /* Wrap pointer around */
  59.     else
  60.         j = r250_index + 103;
  61.  
  62.     new_rand = r250_buffer[r250_index] ^ r250_buffer[j];
  63.     r250_buffer[r250_index] = new_rand;
  64.     if (r250_index >= 249)      /* Increment pointer for next time */
  65.         r250_index = 0;
  66.     else
  67.         r250_index++;
  68.  
  69.     return new_rand;
  70. }
  71. /**** Function:  r250  Description: returns a random double in range 0-1. ****/
  72. double dr250()
  73. {
  74. /*---------------------------------------------------------------------------*/
  75.     register int    j;
  76.     register unsigned int new_rand;
  77. /*---------------------------------------------------------------------------*/
  78.     if (r250_index >= 147)
  79.         j = r250_index - 147;     /* Wrap pointer around */
  80.     else
  81.         j = r250_index + 103;
  82.  
  83.     new_rand = r250_buffer[r250_index] ^ r250_buffer[j];
  84.     r250_buffer[r250_index] = new_rand;
  85.     if (r250_index >= 249)      /* Increment pointer for next time */
  86.         r250_index = 0;
  87.     else
  88.         r250_index++;
  89.     return new_rand / (double)0xffff;   /* Return a number in 0.0 to 1.0 */
  90. }
  91.  
  92.  
  93. [LISTING TWO]
  94.  
  95. /******************************************************************************
  96. *  Module: rtest.c   Description: tests R250 random number generator by 
  97. *  placing data in a set of bins.
  98. ******************************************************************************/
  99.  
  100. #include <stdio.h>
  101. #include <stdlib.h>
  102.  
  103. /**** Constants ****/
  104. #define NMR_RAND    5000
  105. #define MAX_BINS    500
  106.  
  107. /**** Function prototypes *****/
  108. unsigned int r250();
  109. void    r250_init(int seed);
  110.  
  111. /**** Function:     main  ****/
  112. void main(int argc, char *argv[])
  113. {
  114. /*---------------------------------------------------------------------------*/
  115.     int        j, k;
  116.     int        nmr_bins;
  117.     int        seed;
  118.     int        bins[MAX_BINS];
  119.     double    randm;
  120.     double    bin_limit[MAX_BINS];
  121.     double    bin_inc;
  122. /*---------------------------------------------------------------------------*/
  123.     if (argc != 3)
  124.         {
  125.         printf("Usage -- rtest [nmr_bins] [seed]\n");
  126.         exit(1);
  127.         }
  128.     nmr_bins = atoi(argv[1]);
  129.     if (nmr_bins > MAX_BINS)
  130.         {
  131.         printf("Error -- maximum number of bins is %d\n", MAX_BINS);
  132.         exit(1);
  133.         }
  134.     seed = atoi(argv[2]);
  135.     r250_init(seed);
  136.     bin_inc = 1.0 / nmr_bins;
  137.     for (j = 0; j < nmr_bins; j++)
  138.         {
  139.         bins[j] = 0;        // Initialize bins to zero
  140.         bin_limit[j] = (j + 1) * bin_inc;
  141.         }
  142.     bin_limit[nmr_bins-1] = 1.0e7;    // Make sure all others are in last bin
  143.     for (j = 0; j < NMR_RAND; j++)
  144.         {
  145.         randm = r250() / (double)0xffff;
  146.         for (k = 0; k < nmr_bins; k++)
  147.             if (randm < bin_limit[k])
  148.                 {
  149.                 (bins[k])++;
  150.                 break;
  151.                 }
  152.         }
  153.     for (j = 0; j < nmr_bins; j++)
  154.         printf("%d\n", bins[j]);
  155. }
  156.  
  157.