1 #include <u.h> 2 #include <libc.h> 3 4 /* 5 * algorithm by 6 * D. P. Mitchell & J. A. Reeds 7 */ 8 #define LEN 607 9 #define TAP 273 10 #define MASK 0x7fffffffL 11 #define A 48271 12 #define M 2147483647 13 #define Q 44488 14 #define R 3399 15 #define NORM (1.0/(1.0+MASK)) 16 17 static ulong rng_vec[LEN]; 18 static ulong* rng_tap = rng_vec; 19 static ulong* rng_feed = 0; 20 21 void 22 srand(long seed) 23 { 24 long lo, hi, x; 25 int i; 26 27 rng_tap = rng_vec; 28 rng_feed = rng_vec+LEN-TAP; 29 seed = seed%M; 30 if(seed < 0) 31 seed += M; 32 if(seed == 0) 33 seed = 89482311; 34 x = seed; 35 /* 36 * Initialize by x[n+1] = 48271 * x[n] mod (2**31 - 1) 37 */ 38 for(i = -20; i < LEN; i++) { 39 hi = x / Q; 40 lo = x % Q; 41 x = A*lo - R*hi; 42 if(x < 0) 43 x += M; 44 if(i >= 0) 45 rng_vec[i] = x; 46 } 47 } 48 49 long 50 lrand(void) 51 { 52 ulong x; 53 54 rng_tap--; 55 if(rng_tap < rng_vec) { 56 if(rng_feed == 0) { 57 srand(1); 58 rng_tap--; 59 } 60 rng_tap += LEN; 61 } 62 rng_feed--; 63 if(rng_feed < rng_vec) 64 rng_feed += LEN; 65 x = (*rng_feed + *rng_tap) & MASK; 66 *rng_feed = x; 67 return x; 68 } 69 70 int 71 rand(void) 72 { 73 74 return lrand() & 0x7fff; 75 } 76 77 long 78 lnrand(long n) 79 { 80 long slop, v; 81 82 slop = MASK % n; 83 do 84 v = lrand(); 85 while(v <= slop); 86 return v % n; 87 } 88 89 int 90 nrand(int n) 91 { 92 long slop, v; 93 94 slop = MASK % n; 95 do 96 v = lrand(); 97 while(v <= slop); 98 return v % n; 99 } 100 101 double 102 frand(void) 103 { 104 double x; 105 106 do { 107 x = lrand() * NORM; 108 x = (x + lrand()) * NORM; 109 } while(x >= 1); 110 return x; 111 } 112