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
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
|
/* rng.c (pseudo-random number generator) */
/***********************************************************************
* This code is part of GLPK (GNU Linear Programming Kit).
*
* This code is a modified version of the module GB_FLIP, a portable
* pseudo-random number generator. The original version of GB_FLIP is
* a part of The Stanford GraphBase developed by Donald E. Knuth (see
* http://www-cs-staff.stanford.edu/~knuth/sgb.html).
*
* Note that all changes concern only external names, so this modified
* version produces exactly the same results as the original version.
*
* Changes were made by Andrew Makhorin <mao@gnu.org>.
*
* GLPK 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.
*
* GLPK 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 GLPK. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/
#include "env.h"
#include "rng.h"
#if 0
int A[56] = { -1 };
#else
#define A (rand->A)
#endif
/* pseudo-random values */
#if 0
int *fptr = A;
#else
#define fptr (rand->fptr)
#endif
/* the next A value to be exported */
#define mod_diff(x, y) (((x) - (y)) & 0x7FFFFFFF)
/* difference modulo 2^31 */
static int flip_cycle(RNG *rand)
{ /* this is an auxiliary routine to do 55 more steps of the basic
* recurrence, at high speed, and to reset fptr */
int *ii, *jj;
for (ii = &A[1], jj = &A[32]; jj <= &A[55]; ii++, jj++)
*ii = mod_diff(*ii, *jj);
for (jj = &A[1]; ii <= &A[55]; ii++, jj++)
*ii = mod_diff(*ii, *jj);
fptr = &A[54];
return A[55];
}
/***********************************************************************
* NAME
*
* rng_create_rand - create pseudo-random number generator
*
* SYNOPSIS
*
* #include "rng.h"
* RNG *rng_create_rand(void);
*
* DESCRIPTION
*
* The routine rng_create_rand creates and initializes a pseudo-random
* number generator.
*
* RETURNS
*
* The routine returns a pointer to the generator created. */
RNG *rng_create_rand(void)
{ RNG *rand;
int i;
rand = talloc(1, RNG);
A[0] = -1;
for (i = 1; i <= 55; i++) A[i] = 0;
fptr = A;
rng_init_rand(rand, 1);
return rand;
}
/***********************************************************************
* NAME
*
* rng_init_rand - initialize pseudo-random number generator
*
* SYNOPSIS
*
* #include "rng.h"
* void rng_init_rand(RNG *rand, int seed);
*
* DESCRIPTION
*
* The routine rng_init_rand initializes the pseudo-random number
* generator. The parameter seed may be any integer number. Note that
* on creating the generator this routine is called with the parameter
* seed equal to 1. */
void rng_init_rand(RNG *rand, int seed)
{ int i;
int prev = seed, next = 1;
seed = prev = mod_diff(prev, 0);
A[55] = prev;
for (i = 21; i; i = (i + 21) % 55)
{ A[i] = next;
next = mod_diff(prev, next);
if (seed & 1)
seed = 0x40000000 + (seed >> 1);
else
seed >>= 1;
next = mod_diff(next, seed);
prev = A[i];
}
flip_cycle(rand);
flip_cycle(rand);
flip_cycle(rand);
flip_cycle(rand);
flip_cycle(rand);
return;
}
/***********************************************************************
* NAME
*
* rng_next_rand - obtain pseudo-random integer in the range [0, 2^31-1]
*
* SYNOPSIS
*
* #include "rng.h"
* int rng_next_rand(RNG *rand);
*
* RETURNS
*
* The routine rng_next_rand returns a next pseudo-random integer which
* is uniformly distributed between 0 and 2^31-1, inclusive. The period
* length of the generated numbers is 2^85 - 2^30. The low order bits of
* the generated numbers are just as random as the high-order bits. */
int rng_next_rand(RNG *rand)
{ return
*fptr >= 0 ? *fptr-- : flip_cycle(rand);
}
/***********************************************************************
* NAME
*
* rng_unif_rand - obtain pseudo-random integer in the range [0, m-1]
*
* SYNOPSIS
*
* #include "rng.h"
* int rng_unif_rand(RNG *rand, int m);
*
* RETURNS
*
* The routine rng_unif_rand returns a next pseudo-random integer which
* is uniformly distributed between 0 and m-1, inclusive, where m is any
* positive integer less than 2^31. */
#define two_to_the_31 ((unsigned int)0x80000000)
int rng_unif_rand(RNG *rand, int m)
{ unsigned int t = two_to_the_31 - (two_to_the_31 % m);
int r;
xassert(m > 0);
do { r = rng_next_rand(rand); } while (t <= (unsigned int)r);
return r % m;
}
/***********************************************************************
* NAME
*
* rng_delete_rand - delete pseudo-random number generator
*
* SYNOPSIS
*
* #include "rng.h"
* void rng_delete_rand(RNG *rand);
*
* DESCRIPTION
*
* The routine rng_delete_rand frees all the memory allocated to the
* specified pseudo-random number generator. */
void rng_delete_rand(RNG *rand)
{ tfree(rand);
return;
}
/**********************************************************************/
#ifdef GLP_TEST
/* To be sure that this modified version produces the same results as
* the original version, run this validation program. */
int main(void)
{ RNG *rand;
int j;
rand = rng_create_rand();
rng_init_rand(rand, -314159);
if (rng_next_rand(rand) != 119318998)
{ fprintf(stderr, "Failure on the first try!\n");
return -1;
}
for (j = 1; j <= 133; j++) rng_next_rand(rand);
if (rng_unif_rand(rand, 0x55555555) != 748103812)
{ fprintf(stderr, "Failure on the second try!\n");
return -2;
}
fprintf(stderr, "OK, the random-number generator routines seem to"
" work!\n");
rng_delete_rand(rand);
return 0;
}
#endif
/* eof */
|
