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