1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
|
#ifdef HAVE_CONFIG_H
# include <config.h>
#endif
#include "game.hh"
#include <stdio.h>
#include <stdlib.h>
#define M 1000000000
#define M_inv 1.0E-9
#define L 24
#define K 55
#define D 21
/*
Random number generator.
Uses the subtractive series
X[n] = (X[n-55] - X[n-24]) mod M, n >= 55
as suggested in Seminumerical Algorithms, exercise 3.2.2--23.
Knuth gives a FORTRAN implementation in section 3.6 for use as a
portable random number generator.
Michael Coffin
*/
/*
The constant M is machine dependent. It must be even, and
integer arithmetic on the range [-M,+M] must be possible. The
function random() returns values in the half-open interval [0,M).
Minv is just 1/M.
Other possible values for L and K can be chosen from Table 1 of
Section 3.2.2. The constant D should be relatively prime to K and
approximately 0.382*K. (Don't ask me why.)
To get a random number in the real interval [X,Y) the recommended
formula is (X + (frandom() * (Y-X))). For an integral interval,
take the floor.
*/
/* current state of generator */
static u_int32_t X[K]; /* Fibonacci array */
static int cur; /* Current index in array. */
static int curL; /* this is always (cur - L) mod K */
/* seed the generator */
void
zrand_seed(u_int32_t seed)
{
int i;
u_int32_t j, k;
/* Make zrandom() repeatable: Always reset X, cur and curL to their initial
values. -ed */
cur = 0;
curL = K - L;
for (i = 0; i < K; i++)
X[i] = 0;
/* According to Knuth, "This computes a Fibonacci-like sequence; */
/* multiplication of indices by 21 [the constant D] spreads the */
/* values around so as to alleviate initial nonrandomness */
/* problems..." */
seed = seed % M;
X[0] = j = seed;
k = 1;
for (i = 1; i < K; i++) {
int which = (D*i) % K;
X[which] = k;
if (j < k) k -= M;
k = j - k;
j = X[which];
}
/* `warm up' the generator. Three was good enough for Knuth... */
for (i = 0; i < 5 * K; i++)
zrand();
}
/* return next value in the sequence */
u_int32_t
zrand()
{
u_int32_t result;
result = X[cur] < X[curL] ? M : 0;
result += X[cur] - X[curL];
X[cur] = result;
cur++; if (cur == K) cur = 0;
curL++; if (curL == K) curL = 0;
return result;
}
|
