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
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
|
/* rng/knuthran.c
*
* Copyright (C) 2001, 2007 Brian Gough, Carlo Perassi
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or (at
* your option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
/*
* This generator is taken from
*
* Donald E. Knuth
* The Art of Computer Programming
* Volume 2
* Third Edition
* Addison-Wesley
* Section 3.6
*
* The comments are taken from the book
* Our comments are signed
*/
#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>
#define BUFLEN 2009 /* [Brian]: length of the buffer aa[] */
#define KK 100 /* the long lag */
#define LL 37 /* the short lag */
#define MM (1L << 30) /* the modulus */
#define TT 70 /* guaranteed separation between streams */
#define evenize(x) ((x) & (MM - 2)) /* make x even */
#define is_odd(x) ((x) & 1) /* the units bit of x */
#define mod_diff(x, y) (((x) - (y)) & (MM - 1)) /* (x - y) mod MM */
static inline void ran_array (unsigned long int aa[], unsigned int n,
unsigned long int ran_x[]);
static inline unsigned long int ran_get (void *vstate);
static double ran_get_double (void *vstate);
static void ran_set (void *state, unsigned long int s);
typedef struct
{
unsigned int i;
unsigned long int aa[BUFLEN]; /* [Carlo]: I can't pass n to ran_array like
Knuth does */
unsigned long int ran_x[KK]; /* the generator state */
}
ran_state_t;
static inline void
ran_array (unsigned long int aa[], unsigned int n, unsigned long int ran_x[])
{
unsigned int i;
unsigned int j;
for (j = 0; j < KK; j++)
aa[j] = ran_x[j];
for (; j < n; j++)
aa[j] = mod_diff (aa[j - KK], aa[j - LL]);
for (i = 0; i < LL; i++, j++)
ran_x[i] = mod_diff (aa[j - KK], aa[j - LL]);
for (; i < KK; i++, j++)
ran_x[i] = mod_diff (aa[j - KK], ran_x[i - LL]);
}
static inline unsigned long int
ran_get (void *vstate)
{
ran_state_t *state = (ran_state_t *) vstate;
unsigned int i = state->i;
if (i == 0)
{
/* fill buffer with new random numbers */
ran_array (state->aa, BUFLEN, state->ran_x);
}
state->i = (i + 1) % BUFLEN;
return state->aa[i];
}
static double
ran_get_double (void *vstate)
{
ran_state_t *state = (ran_state_t *) vstate;
return ran_get (state) / 1073741824.0; /* [Carlo]: RAND_MAX + 1 */
}
static void
ran_set (void *vstate, unsigned long int s)
{
ran_state_t *state = (ran_state_t *) vstate;
long x[KK + KK - 1]; /* the preparation buffer */
register int j;
register int t;
register long ss = evenize (s + 2);
for (j = 0; j < KK; j++)
{
x[j] = ss; /* bootstrap the buffer */
ss <<= 1;
if (ss >= MM) /* cyclic shift 29 bits */
ss -= MM - 2;
}
for (; j < KK + KK - 1; j++)
x[j] = 0;
x[1]++; /* make x[1] (and only x[1]) odd */
ss = s & (MM - 1);
t = TT - 1;
while (t)
{
for (j = KK - 1; j > 0; j--) /* square */
x[j + j] = x[j];
for (j = KK + KK - 2; j > KK - LL; j -= 2)
x[KK + KK - 1 - j] = evenize (x[j]);
for (j = KK + KK - 2; j >= KK; j--)
if (is_odd (x[j]))
{
x[j - (KK - LL)] = mod_diff (x[j - (KK - LL)], x[j]);
x[j - KK] = mod_diff (x[j - KK], x[j]);
}
if (is_odd (ss))
{ /* multiply by "z" */
for (j = KK; j > 0; j--)
x[j] = x[j - 1];
x[0] = x[KK]; /* shift the buffer cyclically */
if (is_odd (x[KK]))
x[LL] = mod_diff (x[LL], x[KK]);
}
if (ss)
ss >>= 1;
else
t--;
}
state->i = 0;
for (j = 0; j < LL; j++)
state->ran_x[j + KK - LL] = x[j];
for (; j < KK; j++)
state->ran_x[j - LL] = x[j];
return;
}
static const gsl_rng_type ran_type = {
"knuthran", /* name */
0x3fffffffUL, /* RAND_MAX *//* [Carlo]: (2 ^ 30) - 1 */
0, /* RAND_MIN */
sizeof (ran_state_t),
&ran_set,
&ran_get,
&ran_get_double
};
const gsl_rng_type *gsl_rng_knuthran = &ran_type;
|
