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