1 b2cfc4e2 2003-09-30 devnull #include <lib9.h>
4 b2cfc4e2 2003-09-30 devnull * algorithm by
5 b2cfc4e2 2003-09-30 devnull * D. P. Mitchell & J. A. Reeds
8 b2cfc4e2 2003-09-30 devnull #define LEN 607
9 b2cfc4e2 2003-09-30 devnull #define TAP 273
10 b2cfc4e2 2003-09-30 devnull #define MASK 0x7fffffffL
11 b2cfc4e2 2003-09-30 devnull #define A 48271
12 b2cfc4e2 2003-09-30 devnull #define M 2147483647
13 b2cfc4e2 2003-09-30 devnull #define Q 44488
14 b2cfc4e2 2003-09-30 devnull #define R 3399
15 b2cfc4e2 2003-09-30 devnull #define NORM (1.0/(1.0+MASK))
17 b2cfc4e2 2003-09-30 devnull static ulong rng_vec[LEN];
18 b2cfc4e2 2003-09-30 devnull static ulong* rng_tap = rng_vec;
19 b2cfc4e2 2003-09-30 devnull static ulong* rng_feed = 0;
20 b2cfc4e2 2003-09-30 devnull static Lock lk;
22 b2cfc4e2 2003-09-30 devnull static void
23 b2cfc4e2 2003-09-30 devnull isrand(long seed)
25 b2cfc4e2 2003-09-30 devnull long lo, hi, x;
28 b2cfc4e2 2003-09-30 devnull rng_tap = rng_vec;
29 b2cfc4e2 2003-09-30 devnull rng_feed = rng_vec+LEN-TAP;
30 b2cfc4e2 2003-09-30 devnull seed = seed%M;
31 b2cfc4e2 2003-09-30 devnull if(seed < 0)
32 b2cfc4e2 2003-09-30 devnull seed += M;
33 b2cfc4e2 2003-09-30 devnull if(seed == 0)
34 b2cfc4e2 2003-09-30 devnull seed = 89482311;
35 b2cfc4e2 2003-09-30 devnull x = seed;
37 b2cfc4e2 2003-09-30 devnull * Initialize by x[n+1] = 48271 * x[n] mod (2**31 - 1)
39 b2cfc4e2 2003-09-30 devnull for(i = -20; i < LEN; i++) {
40 b2cfc4e2 2003-09-30 devnull hi = x / Q;
41 b2cfc4e2 2003-09-30 devnull lo = x % Q;
42 b2cfc4e2 2003-09-30 devnull x = A*lo - R*hi;
43 b2cfc4e2 2003-09-30 devnull if(x < 0)
45 b2cfc4e2 2003-09-30 devnull if(i >= 0)
46 b2cfc4e2 2003-09-30 devnull rng_vec[i] = x;
51 b2cfc4e2 2003-09-30 devnull srand(long seed)
53 b2cfc4e2 2003-09-30 devnull lock(&lk);
54 b2cfc4e2 2003-09-30 devnull isrand(seed);
55 b2cfc4e2 2003-09-30 devnull unlock(&lk);
59 b2cfc4e2 2003-09-30 devnull lrand(void)
63 b2cfc4e2 2003-09-30 devnull lock(&lk);
65 b2cfc4e2 2003-09-30 devnull rng_tap--;
66 b2cfc4e2 2003-09-30 devnull if(rng_tap < rng_vec) {
67 b2cfc4e2 2003-09-30 devnull if(rng_feed == 0) {
68 b2cfc4e2 2003-09-30 devnull isrand(1);
69 b2cfc4e2 2003-09-30 devnull rng_tap--;
71 b2cfc4e2 2003-09-30 devnull rng_tap += LEN;
73 b2cfc4e2 2003-09-30 devnull rng_feed--;
74 b2cfc4e2 2003-09-30 devnull if(rng_feed < rng_vec)
75 b2cfc4e2 2003-09-30 devnull rng_feed += LEN;
76 b2cfc4e2 2003-09-30 devnull x = (*rng_feed + *rng_tap) & MASK;
77 b2cfc4e2 2003-09-30 devnull *rng_feed = x;
79 b2cfc4e2 2003-09-30 devnull unlock(&lk);
81 b2cfc4e2 2003-09-30 devnull return x;
85 b2cfc4e2 2003-09-30 devnull rand(void)
87 b2cfc4e2 2003-09-30 devnull return lrand() & 0x7fff;