BGI/soap_src/soap_builder/r250.c

   1 /* r250.c       the r250 uniform random number algorithm
   2
   3                 Kirkpatrick, S., and E. Stoll, 1981; "A Very Fast
   4                 Shift-Register Sequence Random Number Generator",
   5                 Journal of Computational Physics, V.40
   6
   7                 also:
   8
   9                 see W.L. Maier, DDJ May 1991
  10
  11
  12
  13 */
  14
  15 #include <limits.h>
  16 #include "r250.h"
  17
  18 // static functions
  19 static unsigned int randlcg();
  20
  21
  22 #define MSB          0x80000000L
  23 #define ALL_BITS     0xffffffffL
  24 #define HALF_RANGE   0x40000000L
  25 #define STEP         7
  26 #define BITS         32
  27
  28 static unsigned int r250_buffer[ 250 ];
  29 static int r250_index;
  30
  31 static unsigned int randlcg(int sd)       /* returns a random unsigned integer */
  32 {
  33                 static long quotient1  = LONG_MAX / 16807L;
  34                 static long remainder1 = LONG_MAX % 16807L;
  35
  36         if ( sd <= quotient1 )
  37                 sd = (sd * 16807L) % LONG_MAX;
  38         else
  39         {
  40                 int high_part = sd / quotient1;
  41                 int low_part  = sd % quotient1;
  42
  43                 int test = 16807L * low_part - remainder1 * high_part;
  44
  45                 if ( test > 0 )
  46                         sd = test;
  47                 else
  48                         sd = test + LONG_MAX;
  49
  50         }
  51
  52         return sd;
  53 }
  54
  55 void r250_init(int sd)
  56 {
  57         int j, k;
  58         unsigned int mask, msb;
  59
  60         r250_index = 0;
  61         for (j = 0; j < 250; j++)      /* fill r250 buffer with BITS-1 bit values */
  62                 sd = r250_buffer[j] = randlcg(sd);
  63
  64
  65         for (j = 0; j < 250; j++)       /* set some MSBs to 1 */
  66                 if ( (sd=randlcg(sd)) > HALF_RANGE )
  67                         r250_buffer[j] |= MSB;
  68
  69
  70         msb = MSB;              /* turn on diagonal bit */
  71         mask = ALL_BITS;        /* turn off the leftmost bits */
  72
  73         for (j = 0; j < BITS; j++)
  74         {
  75                 k = STEP * j + 3;       /* select a word to operate on */
  76                 r250_buffer[k] &= mask; /* turn off bits left of the diagonal */
  77                 r250_buffer[k] |= msb;  /* turn on the diagonal bit */
  78                 mask >>= 1;
  79                 msb  >>= 1;
  80         }
  81
  82 }
  83
  84 unsigned int r250()             /* returns a random unsigned integer */
  85 {
  86         register int    j;
  87         register unsigned int new_rand;
  88
  89         if ( r250_index >= 147 )
  90                 j = r250_index - 147;   /* wrap pointer around */
  91         else
  92                 j = r250_index + 103;
  93
  94         new_rand = r250_buffer[ r250_index ] ^ r250_buffer[ j ];
  95         r250_buffer[ r250_index ] = new_rand;
  96
  97         if ( r250_index >= 249 )        /* increment pointer for next time */
  98                 r250_index = 0;
  99         else
 100                 r250_index++;
 101
 102         return new_rand;
 103
 104 }
 105
 106
 107 double dr250()          /* returns a random double in range 0..1 */
 108 {
 109         register int    j;
 110         register unsigned int new_rand;
 111
 112         if ( r250_index >= 147 )
 113                 j = r250_index - 147;   /* wrap pointer around */
 114         else
 115                 j = r250_index + 103;
 116
 117         new_rand = r250_buffer[ r250_index ] ^ r250_buffer[ j ];
 118         r250_buffer[ r250_index ] = new_rand;
 119
 120         if ( r250_index >= 249 )        /* increment pointer for next time */
 121                 r250_index = 0;
 122         else
 123                 r250_index++;
 124
 125         return (double)new_rand / ALL_BITS;
 126
 127 }
 128