1 b2cfc4e2 2003-09-30 devnull #include	<lib9.h>
 2 b2cfc4e2 2003-09-30 devnull 
 3 b2cfc4e2 2003-09-30 devnull /*
 4 b2cfc4e2 2003-09-30 devnull  *	algorithm by
 5 b2cfc4e2 2003-09-30 devnull  *	D. P. Mitchell & J. A. Reeds
 6 b2cfc4e2 2003-09-30 devnull  */
 7 b2cfc4e2 2003-09-30 devnull 
 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))
16 b2cfc4e2 2003-09-30 devnull 
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;
21 b2cfc4e2 2003-09-30 devnull 
22 b2cfc4e2 2003-09-30 devnull static void
23 b2cfc4e2 2003-09-30 devnull isrand(long seed)
24 b2cfc4e2 2003-09-30 devnull {
25 b2cfc4e2 2003-09-30 devnull 	long lo, hi, x;
26 b2cfc4e2 2003-09-30 devnull 	int i;
27 b2cfc4e2 2003-09-30 devnull 
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;
36 b2cfc4e2 2003-09-30 devnull 	/*
37 b2cfc4e2 2003-09-30 devnull 	 *	Initialize by x[n+1] = 48271 * x[n] mod (2**31 - 1)
38 b2cfc4e2 2003-09-30 devnull 	 */
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)
44 b2cfc4e2 2003-09-30 devnull 			x += M;
45 b2cfc4e2 2003-09-30 devnull 		if(i >= 0)
46 b2cfc4e2 2003-09-30 devnull 			rng_vec[i] = x;
47 b2cfc4e2 2003-09-30 devnull 	}
48 b2cfc4e2 2003-09-30 devnull }
49 b2cfc4e2 2003-09-30 devnull 
50 b2cfc4e2 2003-09-30 devnull void
51 b2cfc4e2 2003-09-30 devnull srand(long seed)
52 b2cfc4e2 2003-09-30 devnull {
53 b2cfc4e2 2003-09-30 devnull 	lock(&lk);
54 b2cfc4e2 2003-09-30 devnull 	isrand(seed);
55 b2cfc4e2 2003-09-30 devnull 	unlock(&lk);
56 b2cfc4e2 2003-09-30 devnull }
57 b2cfc4e2 2003-09-30 devnull 
58 b2cfc4e2 2003-09-30 devnull long
59 b2cfc4e2 2003-09-30 devnull lrand(void)
60 b2cfc4e2 2003-09-30 devnull {
61 b2cfc4e2 2003-09-30 devnull 	ulong x;
62 b2cfc4e2 2003-09-30 devnull 
63 b2cfc4e2 2003-09-30 devnull 	lock(&lk);
64 b2cfc4e2 2003-09-30 devnull 
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--;
70 b2cfc4e2 2003-09-30 devnull 		}
71 b2cfc4e2 2003-09-30 devnull 		rng_tap += LEN;
72 b2cfc4e2 2003-09-30 devnull 	}
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;
78 b2cfc4e2 2003-09-30 devnull 
79 b2cfc4e2 2003-09-30 devnull 	unlock(&lk);
80 b2cfc4e2 2003-09-30 devnull 
81 b2cfc4e2 2003-09-30 devnull 	return x;
82 b2cfc4e2 2003-09-30 devnull }
83 b2cfc4e2 2003-09-30 devnull 
84 b2cfc4e2 2003-09-30 devnull int
85 b2cfc4e2 2003-09-30 devnull rand(void)
86 b2cfc4e2 2003-09-30 devnull {
87 b2cfc4e2 2003-09-30 devnull 	return lrand() & 0x7fff;
88 b2cfc4e2 2003-09-30 devnull }
89 b2cfc4e2 2003-09-30 devnull