13e12c5d1SDavid du Colombier #include <stdlib.h> 23e12c5d1SDavid du Colombier 33e12c5d1SDavid du Colombier /* 43e12c5d1SDavid du Colombier * algorithm by 53e12c5d1SDavid du Colombier * D. P. Mitchell & J. A. Reeds 63e12c5d1SDavid du Colombier */ 73e12c5d1SDavid du Colombier #define LEN 607 83e12c5d1SDavid du Colombier #define TAP 273 93e12c5d1SDavid du Colombier #define MASK 0x7fffffffL 103e12c5d1SDavid du Colombier #define A 48271 113e12c5d1SDavid du Colombier #define M 2147483647 123e12c5d1SDavid du Colombier #define Q 44488 133e12c5d1SDavid du Colombier #define R 3399 143e12c5d1SDavid du Colombier 153e12c5d1SDavid du Colombier typedef unsigned long ulong; 163e12c5d1SDavid du Colombier 173e12c5d1SDavid du Colombier static ulong rng_vec[LEN]; 183e12c5d1SDavid du Colombier static ulong* rng_tap = rng_vec; 193e12c5d1SDavid du Colombier static ulong* rng_feed = 0; 203e12c5d1SDavid du Colombier 213e12c5d1SDavid du Colombier void srand(unsigned int seed)223e12c5d1SDavid du Colombiersrand(unsigned int seed) 233e12c5d1SDavid du Colombier { 243e12c5d1SDavid du Colombier long lo, hi, x; 253e12c5d1SDavid du Colombier int i; 263e12c5d1SDavid du Colombier 273e12c5d1SDavid du Colombier rng_tap = rng_vec; 283e12c5d1SDavid du Colombier rng_feed = rng_vec+LEN-TAP; 293e12c5d1SDavid du Colombier seed = seed%M; 30*781103c4SDavid du Colombier if(0 && seed < 0) /* seed is unsigned */ 313e12c5d1SDavid du Colombier seed += M; 323e12c5d1SDavid du Colombier if(seed == 0) 333e12c5d1SDavid du Colombier seed = 89482311; 343e12c5d1SDavid du Colombier x = seed; 353e12c5d1SDavid du Colombier /* 363e12c5d1SDavid du Colombier * Initialize by x[n+1] = 48271 * x[n] mod (2**31 - 1) 373e12c5d1SDavid du Colombier */ 383e12c5d1SDavid du Colombier for(i = -20; i < LEN; i++) { 393e12c5d1SDavid du Colombier hi = x / Q; 403e12c5d1SDavid du Colombier lo = x % Q; 413e12c5d1SDavid du Colombier x = A*lo - R*hi; 423e12c5d1SDavid du Colombier if(x < 0) 433e12c5d1SDavid du Colombier x += M; 443e12c5d1SDavid du Colombier if(i >= 0) 453e12c5d1SDavid du Colombier rng_vec[i] = x; 463e12c5d1SDavid du Colombier } 473e12c5d1SDavid du Colombier } 483e12c5d1SDavid du Colombier 493e12c5d1SDavid du Colombier static long lrand(void)503e12c5d1SDavid du Colombierlrand(void) 513e12c5d1SDavid du Colombier { 523e12c5d1SDavid du Colombier ulong x; 533e12c5d1SDavid du Colombier 543e12c5d1SDavid du Colombier rng_tap--; 553e12c5d1SDavid du Colombier if(rng_tap < rng_vec) { 563e12c5d1SDavid du Colombier if(rng_feed == 0) { 573e12c5d1SDavid du Colombier srand(1); 583e12c5d1SDavid du Colombier rng_tap--; 593e12c5d1SDavid du Colombier } 603e12c5d1SDavid du Colombier rng_tap += LEN; 613e12c5d1SDavid du Colombier } 623e12c5d1SDavid du Colombier rng_feed--; 633e12c5d1SDavid du Colombier if(rng_feed < rng_vec) 643e12c5d1SDavid du Colombier rng_feed += LEN; 653e12c5d1SDavid du Colombier x = (*rng_feed + *rng_tap) & MASK; 663e12c5d1SDavid du Colombier *rng_feed = x; 673e12c5d1SDavid du Colombier return x; 683e12c5d1SDavid du Colombier } 693e12c5d1SDavid du Colombier 703e12c5d1SDavid du Colombier int rand(void)713e12c5d1SDavid du Colombierrand(void) 723e12c5d1SDavid du Colombier { 733e12c5d1SDavid du Colombier 743e12c5d1SDavid du Colombier return lrand() & 0x7fff; 753e12c5d1SDavid du Colombier } 76