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