sys/src/cmd/cwfs/lrand.c

   1 #include "all.h"
   2
   3 /*
   4  *      algorithm by
   5  *      D. P. Mitchell & J. A. Reeds
   6  */
   7 enum {
   8         LEN     = 607,
   9         TAP     = 273,
  10         MASK    = 0x7fffffffL,
  11         A       = 48271,
  12         M       = 2147483647,
  13         Q       = 44488,
  14         R       = 3399,
  15 };
  16
  17 #define NORM    (1.0/(1.0+MASK))
  18
  19 static  ulong   rng_vec[LEN];
  20 static  ulong*  rng_tap = rng_vec;
  21 static  ulong*  rng_feed = 0;
  22 static  Lock    lk;
  23
  24 static void
  25 isrand(long seed)
  26 {
  27         long lo, hi, x;
  28         int i;
  29
  30         rng_tap = rng_vec;
  31         rng_feed = rng_vec+LEN-TAP;
  32         seed %= M;
  33         if(seed < 0)
  34                 seed += M;
  35         if(seed == 0)
  36                 seed = 89482311;
  37         x = seed;
  38         /*
  39          *      Initialize by x[n+1] = 48271 * x[n] mod (2**31 - 1)
  40          */
  41         for(i = -20; i < LEN; i++) {
  42                 hi = x / Q;
  43                 lo = x % Q;
  44                 x = A*lo - R*hi;
  45                 if(x < 0)
  46                         x += M;
  47                 if(i >= 0)
  48                         rng_vec[i] = x;
  49         }
  50 }
  51
  52 void
  53 srand(long seed)
  54 {
  55         lock(&lk);
  56         isrand(seed);
  57         unlock(&lk);
  58 }
  59
  60 long
  61 lrand(void)
  62 {
  63         ulong x;
  64
  65         lock(&lk);
  66
  67         rng_tap--;
  68         if(rng_tap < rng_vec) {
  69                 if(rng_feed == 0) {
  70                         isrand(1);
  71                         rng_tap--;
  72                 }
  73                 rng_tap += LEN;
  74         }
  75         rng_feed--;
  76         if(rng_feed < rng_vec)
  77                 rng_feed += LEN;
  78         x = (*rng_feed + *rng_tap) & MASK;
  79         *rng_feed = x;
  80
  81         unlock(&lk);
  82         return x;
  83 }