1 /* rng.c (pseudo-random number generator) */
2
3 /***********************************************************************
4 * This code is part of GLPK (GNU Linear Programming Kit).
5 *
6 * This code is a modified version of the module GB_FLIP, a portable
7 * pseudo-random number generator. The original version of GB_FLIP is
8 * a part of The Stanford GraphBase developed by Donald E. Knuth (see
9 * http://www-cs-staff.stanford.edu/~knuth/sgb.html).
10 *
11 * Note that all changes concern only external names, so this modified
12 * version produces exactly the same results as the original version.
13 *
14 * Changes were made by Andrew Makhorin <mao@gnu.org>.
15 *
16 * GLPK is free software: you can redistribute it and/or modify it
17 * under the terms of the GNU General Public License as published by
18 * the Free Software Foundation, either version 3 of the License, or
19 * (at your option) any later version.
20 *
21 * GLPK is distributed in the hope that it will be useful, but WITHOUT
22 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
23 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
24 * License for more details.
25 *
26 * You should have received a copy of the GNU General Public License
27 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
28 ***********************************************************************/
29
30 #include "env.h"
31 #include "rng.h"
32
33 #if 0
34 int A[56] = { -1 };
35 #else
36 #define A (rand->A)
37 #endif
38 /* pseudo-random values */
39
40 #if 0
41 int *fptr = A;
42 #else
43 #define fptr (rand->fptr)
44 #endif
45 /* the next A value to be exported */
46
47 #define mod_diff(x, y) (((x) - (y)) & 0x7FFFFFFF)
48 /* difference modulo 2^31 */
49
flip_cycle(RNG * rand)50 static int flip_cycle(RNG *rand)
51 { /* this is an auxiliary routine to do 55 more steps of the basic
52 * recurrence, at high speed, and to reset fptr */
53 int *ii, *jj;
54 for (ii = &A[1], jj = &A[32]; jj <= &A[55]; ii++, jj++)
55 *ii = mod_diff(*ii, *jj);
56 for (jj = &A[1]; ii <= &A[55]; ii++, jj++)
57 *ii = mod_diff(*ii, *jj);
58 fptr = &A[54];
59 return A[55];
60 }
61
62 /***********************************************************************
63 * NAME
64 *
65 * rng_create_rand - create pseudo-random number generator
66 *
67 * SYNOPSIS
68 *
69 * #include "rng.h"
70 * RNG *rng_create_rand(void);
71 *
72 * DESCRIPTION
73 *
74 * The routine rng_create_rand creates and initializes a pseudo-random
75 * number generator.
76 *
77 * RETURNS
78 *
79 * The routine returns a pointer to the generator created. */
80
rng_create_rand(void)81 RNG *rng_create_rand(void)
82 { RNG *rand;
83 int i;
84 rand = talloc(1, RNG);
85 A[0] = -1;
86 for (i = 1; i <= 55; i++) A[i] = 0;
87 fptr = A;
88 rng_init_rand(rand, 1);
89 return rand;
90 }
91
92 /***********************************************************************
93 * NAME
94 *
95 * rng_init_rand - initialize pseudo-random number generator
96 *
97 * SYNOPSIS
98 *
99 * #include "rng.h"
100 * void rng_init_rand(RNG *rand, int seed);
101 *
102 * DESCRIPTION
103 *
104 * The routine rng_init_rand initializes the pseudo-random number
105 * generator. The parameter seed may be any integer number. Note that
106 * on creating the generator this routine is called with the parameter
107 * seed equal to 1. */
108
rng_init_rand(RNG * rand,int seed)109 void rng_init_rand(RNG *rand, int seed)
110 { int i;
111 int prev = seed, next = 1;
112 seed = prev = mod_diff(prev, 0);
113 A[55] = prev;
114 for (i = 21; i; i = (i + 21) % 55)
115 { A[i] = next;
116 next = mod_diff(prev, next);
117 if (seed & 1)
118 seed = 0x40000000 + (seed >> 1);
119 else
120 seed >>= 1;
121 next = mod_diff(next, seed);
122 prev = A[i];
123 }
124 flip_cycle(rand);
125 flip_cycle(rand);
126 flip_cycle(rand);
127 flip_cycle(rand);
128 flip_cycle(rand);
129 return;
130 }
131
132 /***********************************************************************
133 * NAME
134 *
135 * rng_next_rand - obtain pseudo-random integer in the range [0, 2^31-1]
136 *
137 * SYNOPSIS
138 *
139 * #include "rng.h"
140 * int rng_next_rand(RNG *rand);
141 *
142 * RETURNS
143 *
144 * The routine rng_next_rand returns a next pseudo-random integer which
145 * is uniformly distributed between 0 and 2^31-1, inclusive. The period
146 * length of the generated numbers is 2^85 - 2^30. The low order bits of
147 * the generated numbers are just as random as the high-order bits. */
148
rng_next_rand(RNG * rand)149 int rng_next_rand(RNG *rand)
150 { return
151 *fptr >= 0 ? *fptr-- : flip_cycle(rand);
152 }
153
154 /***********************************************************************
155 * NAME
156 *
157 * rng_unif_rand - obtain pseudo-random integer in the range [0, m-1]
158 *
159 * SYNOPSIS
160 *
161 * #include "rng.h"
162 * int rng_unif_rand(RNG *rand, int m);
163 *
164 * RETURNS
165 *
166 * The routine rng_unif_rand returns a next pseudo-random integer which
167 * is uniformly distributed between 0 and m-1, inclusive, where m is any
168 * positive integer less than 2^31. */
169
170 #define two_to_the_31 ((unsigned int)0x80000000)
171
rng_unif_rand(RNG * rand,int m)172 int rng_unif_rand(RNG *rand, int m)
173 { unsigned int t = two_to_the_31 - (two_to_the_31 % m);
174 int r;
175 xassert(m > 0);
176 do { r = rng_next_rand(rand); } while (t <= (unsigned int)r);
177 return r % m;
178 }
179
180 /***********************************************************************
181 * NAME
182 *
183 * rng_delete_rand - delete pseudo-random number generator
184 *
185 * SYNOPSIS
186 *
187 * #include "rng.h"
188 * void rng_delete_rand(RNG *rand);
189 *
190 * DESCRIPTION
191 *
192 * The routine rng_delete_rand frees all the memory allocated to the
193 * specified pseudo-random number generator. */
194
rng_delete_rand(RNG * rand)195 void rng_delete_rand(RNG *rand)
196 { tfree(rand);
197 return;
198 }
199
200 /**********************************************************************/
201
202 #ifdef GLP_TEST
203 /* To be sure that this modified version produces the same results as
204 * the original version, run this validation program. */
205
main(void)206 int main(void)
207 { RNG *rand;
208 int j;
209 rand = rng_create_rand();
210 rng_init_rand(rand, -314159);
211 if (rng_next_rand(rand) != 119318998)
212 { fprintf(stderr, "Failure on the first try!\n");
213 return -1;
214 }
215 for (j = 1; j <= 133; j++) rng_next_rand(rand);
216 if (rng_unif_rand(rand, 0x55555555) != 748103812)
217 { fprintf(stderr, "Failure on the second try!\n");
218 return -2;
219 }
220 fprintf(stderr, "OK, the random-number generator routines seem to"
221 " work!\n");
222 rng_delete_rand(rand);
223 return 0;
224 }
225 #endif
226
227 /* eof */
228