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