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