Subversion Repositories planix.SVN

Rev

Rev 2 | Details | Compare with Previous | Last modification | View Log | RSS feed

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