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 isrand(long seed)25isrand(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 srand(long seed)53srand(long seed) 54 { 55 lock(&lk); 56 isrand(seed); 57 unlock(&lk); 58 } 59 60 long lrand(void)61lrand(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 } 84