Bug Summary

File:gromacs/random/random.c
Location:line 264, column 5
Description:Value stored to 'mti' is never read

Annotated Source Code

1/*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2008, The GROMACS development team.
6 * Copyright (c) 2012,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
10 *
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
15 *
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
20 *
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 *
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
33 *
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
36 */
37#ifdef HAVE_CONFIG_H1
38#include <config.h>
39#endif
40
41#include "random.h"
42
43#include <stdio.h>
44#include <stdlib.h>
45#ifdef HAVE_UNISTD_H
46#include <unistd.h>
47#endif
48#include <time.h>
49#include <math.h>
50#ifdef GMX_NATIVE_WINDOWS
51#include <process.h>
52#endif
53
54#include "external/Random123-1.08/include/Random123/threefry.h"
55
56#include "gromacs/math/utilities.h"
57#include "random_gausstable.h"
58
59#define RNG_N624 624
60#define RNG_M397 397
61#define RNG_MATRIX_A0x9908b0dfUL 0x9908b0dfUL /* constant vector a */
62#define RNG_UPPER_MASK0x80000000UL 0x80000000UL /* most significant w-r bits */
63#define RNG_LOWER_MASK0x7fffffffUL 0x7fffffffUL /* least significant r bits */
64
65/* Note that if you change the size of the Gaussian table you will
66 * also have to generate new initialization data for the table in
67 * gmx_random_gausstable.h
68 *
69 * A routine print_gaussian_table() is in contrib/random.c
70 * for convenience - use it if you need a different size of the table.
71 */
72#define GAUSS_TABLE14 14 /* the size of the gauss table is 2^GAUSS_TABLE */
73#define GAUSS_MASK((1 << 14) - 1) ((1 << GAUSS_TABLE14) - 1)
74
75
76struct gmx_rng {
77 unsigned int mt[RNG_N624];
78 int mti;
79 int has_saved;
80 double gauss_saved;
81};
82
83
84
85int
86gmx_rng_n(void)
87{
88 return RNG_N624;
89}
90
91
92gmx_rng_t
93gmx_rng_init(unsigned int seed)
94{
95 struct gmx_rng *rng;
96
97 if ((rng = (struct gmx_rng *)malloc(sizeof(struct gmx_rng))) == NULL((void*)0))
98 {
99 return NULL((void*)0);
100 }
101
102 rng->has_saved = 0; /* no saved gaussian number yet */
103
104 rng->mt[0] = seed & 0xffffffffUL;
105 for (rng->mti = 1; rng->mti < RNG_N624; rng->mti++)
106 {
107 rng->mt[rng->mti] =
108 (1812433253UL * (rng->mt[rng->mti-1] ^
109 (rng->mt[rng->mti-1] >> 30)) + rng->mti);
110 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
111 /* In the previous versions, MSBs of the seed affect */
112 /* only MSBs of the array mt[]. */
113 /* 2002/01/09 modified by Makoto Matsumoto */
114 rng->mt[rng->mti] &= 0xffffffffUL;
115 /* for >32 bit machines */
116 }
117 return rng;
118}
119
120gmx_rng_t
121gmx_rng_init_array(unsigned int seed[], int seed_length)
122{
123 int i, j, k;
124 gmx_rng_t rng;
125
126 if ((rng = gmx_rng_init(19650218UL)) == NULL((void*)0))
127 {
128 return NULL((void*)0);
129 }
130
131 i = 1; j = 0;
132 k = (RNG_N624 > seed_length ? RNG_N624 : seed_length);
133 for (; k; k--)
134 {
135 rng->mt[i] = (rng->mt[i] ^ ((rng->mt[i-1] ^
136 (rng->mt[i-1] >> 30)) * 1664525UL))
137 + seed[j] + j; /* non linear */
138 rng->mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
139 i++; j++;
140 if (i >= RNG_N624)
141 {
142 rng->mt[0] = rng->mt[RNG_N624-1]; i = 1;
143 }
144 if (j >= seed_length)
145 {
146 j = 0;
147 }
148 }
149 for (k = RNG_N624-1; k; k--)
150 {
151 rng->mt[i] = (rng->mt[i] ^ ((rng->mt[i-1] ^
152 (rng->mt[i-1] >> 30)) *
153 1566083941UL))
154 - i; /* non linear */
155 rng->mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
156 i++;
157 if (i >= RNG_N624)
158 {
159 rng->mt[0] = rng->mt[RNG_N624-1]; i = 1;
160 }
161 }
162
163 rng->mt[0] = 0x80000000UL;
164 /* MSB is 1; assuring non-zero initial array */
165 return rng;
166}
167
168
169void
170gmx_rng_destroy(gmx_rng_t rng)
171{
172 if (rng)
173 {
174 free(rng);
175 }
176}
177
178
179void
180gmx_rng_get_state(gmx_rng_t rng, unsigned int *mt, int *mti)
181{
182 int i;
183
184 for (i = 0; i < RNG_N624; i++)
185 {
186 mt[i] = rng->mt[i];
187 }
188 *mti = rng->mti;
189}
190
191
192void
193gmx_rng_set_state(gmx_rng_t rng, unsigned int *mt, int mti)
194{
195 int i;
196
197 for (i = 0; i < RNG_N624; i++)
198 {
199 rng->mt[i] = mt[i];
200 }
201 rng->mti = mti;
202}
203
204
205unsigned int
206gmx_rng_make_seed(void)
207{
208 FILE *fp;
209 unsigned int data;
210 int ret;
211 long my_pid;
212
213#ifdef HAVE_UNISTD_H
214 /* We never want Gromacs execution to halt 10-20 seconds while
215 * waiting for enough entropy in the random number generator.
216 * For this reason we should NOT use /dev/random, which will
217 * block in cases like that. That will cause all sorts of
218 * Gromacs programs to block ~20 seconds while waiting for a
219 * super-random-number to generate cool quotes. Apart from the
220 * minor irritation, it is really bad behavior of a program
221 * to abuse the system random numbers like that - other programs
222 * need them too.
223 * For this reason, we ONLY try to get random numbers from
224 * the pseudo-random stream /dev/urandom, and use other means
225 * if it is not present (in this case fopen() returns NULL).
226 */
227 fp = fopen("/dev/urandom", "rb");
228#else
229 fp = NULL((void*)0);
230#endif
231 if (fp != NULL((void*)0))
232 {
233 ret = fread(&data, sizeof(unsigned int), 1, fp);
234 fclose(fp);
235 }
236 else
237 {
238 /* No random device available, use time-of-day and process id */
239#ifdef GMX_NATIVE_WINDOWS
240 my_pid = (long)_getpid();
241#else
242 my_pid = (long)getpid();
243#endif
244 data = (unsigned int)(((long)time(NULL((void*)0))+my_pid) % (long)1000000);
245 }
246 return data;
247}
248
249
250/* The random number state contains RNG_N entries that are returned one by
251 * one as random numbers. When we run out of them, this routine is called to
252 * regenerate RNG_N new entries.
253 */
254static void
255gmx_rng_update(gmx_rng_t rng)
256{
257 unsigned int lastx, x1, x2, y, *mt;
258 int mti, k;
259 const unsigned int mag01[2] = {0x0UL, RNG_MATRIX_A0x9908b0dfUL};
260 /* mag01[x] = x * MATRIX_A for x=0,1 */
261
262 /* update random numbers */
263 mt = rng->mt; /* pointer to array - avoid repeated dereferencing */
264 mti = rng->mti;
Value stored to 'mti' is never read
265
266 x1 = mt[0];
267 for (k = 0; k < RNG_N624-RNG_M397-3; k += 4)
268 {
269 x2 = mt[k+1];
270 y = (x1 & RNG_UPPER_MASK0x80000000UL) | (x2 & RNG_LOWER_MASK0x7fffffffUL);
271 mt[k] = mt[k+RNG_M397] ^ (y >> 1) ^ mag01[y & 0x1UL];
272 x1 = mt[k+2];
273 y = (x2 & RNG_UPPER_MASK0x80000000UL) | (x1 & RNG_LOWER_MASK0x7fffffffUL);
274 mt[k+1] = mt[k+RNG_M397+1] ^ (y >> 1) ^ mag01[y & 0x1UL];
275 x2 = mt[k+3];
276 y = (x1 & RNG_UPPER_MASK0x80000000UL) | (x2 & RNG_LOWER_MASK0x7fffffffUL);
277 mt[k+2] = mt[k+RNG_M397+2] ^ (y >> 1) ^ mag01[y & 0x1UL];
278 x1 = mt[k+4];
279 y = (x2 & RNG_UPPER_MASK0x80000000UL) | (x1 & RNG_LOWER_MASK0x7fffffffUL);
280 mt[k+3] = mt[k+RNG_M397+3] ^ (y >> 1) ^ mag01[y & 0x1UL];
281 }
282 x2 = mt[k+1];
283 y = (x1 & RNG_UPPER_MASK0x80000000UL) | (x2 & RNG_LOWER_MASK0x7fffffffUL);
284 mt[k] = mt[k+RNG_M397] ^ (y >> 1) ^ mag01[y & 0x1UL];
285 k++;
286 x1 = mt[k+1];
287 y = (x2 & RNG_UPPER_MASK0x80000000UL) | (x1 & RNG_LOWER_MASK0x7fffffffUL);
288 mt[k] = mt[k+RNG_M397] ^ (y >> 1) ^ mag01[y & 0x1UL];
289 k++;
290 x2 = mt[k+1];
291 y = (x1 & RNG_UPPER_MASK0x80000000UL) | (x2 & RNG_LOWER_MASK0x7fffffffUL);
292 mt[k] = mt[k+RNG_M397] ^ (y >> 1) ^ mag01[y & 0x1UL];
293 k++;
294 for (; k < RNG_N624-1; k += 4)
295 {
296 x1 = mt[k+1];
297 y = (x2 & RNG_UPPER_MASK0x80000000UL) | (x1 & RNG_LOWER_MASK0x7fffffffUL);
298 mt[k] = mt[k+(RNG_M397-RNG_N624)] ^ (y >> 1) ^ mag01[y & 0x1UL];
299 x2 = mt[k+2];
300 y = (x1 & RNG_UPPER_MASK0x80000000UL) | (x2 & RNG_LOWER_MASK0x7fffffffUL);
301 mt[k+1] = mt[k+(RNG_M397-RNG_N624)+1] ^ (y >> 1) ^ mag01[y & 0x1UL];
302 x1 = mt[k+3];
303 y = (x2 & RNG_UPPER_MASK0x80000000UL) | (x1 & RNG_LOWER_MASK0x7fffffffUL);
304 mt[k+2] = mt[k+(RNG_M397-RNG_N624)+2] ^ (y >> 1) ^ mag01[y & 0x1UL];
305 x2 = mt[k+4];
306 y = (x1 & RNG_UPPER_MASK0x80000000UL) | (x2 & RNG_LOWER_MASK0x7fffffffUL);
307 mt[k+3] = mt[k+(RNG_M397-RNG_N624)+3] ^ (y >> 1) ^ mag01[y & 0x1UL];
308 }
309 y = (x2 & RNG_UPPER_MASK0x80000000UL) | (mt[0] & RNG_LOWER_MASK0x7fffffffUL);
310 mt[RNG_N624-1] = mt[RNG_M397-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
311
312 rng->mti = 0;
313}
314
315
316real
317gmx_rng_gaussian_real(gmx_rng_t rng)
318{
319 real x, y, r;
320
321 if (rng->has_saved)
322 {
323 rng->has_saved = 0;
324 return rng->gauss_saved;
325 }
326 else
327 {
328 do
329 {
330 x = 2.0*gmx_rng_uniform_real(rng)-1.0;
331 y = 2.0*gmx_rng_uniform_real(rng)-1.0;
332 r = x*x+y*y;
333 }
334 while (r > 1.0 || r == 0.0);
335
336 r = sqrt(-2.0*log(r)/r);
337 rng->gauss_saved = y*r; /* save second random number */
338 rng->has_saved = 1;
339 return x*r; /* return first random number */
340 }
341}
342
343
344
345
346/* Return a random unsigned integer, i.e. 0..4294967295
347 * Provided in header file for performace reasons.
348 * Unfortunately this function cannot be inlined, since
349 * it needs to refer the internal-linkage gmx_rng_update().
350 */
351unsigned int
352gmx_rng_uniform_uint32(gmx_rng_t rng)
353{
354 unsigned int y;
355
356 if (rng->mti == RNG_N624)
357 {
358 gmx_rng_update(rng);
359 }
360 y = rng->mt[rng->mti++];
361
362 y ^= (y >> 11);
363 y ^= (y << 7) & 0x9d2c5680UL;
364 y ^= (y << 15) & 0xefc60000UL;
365 y ^= (y >> 18);
366
367 return y;
368}
369
370
371
372
373
374/* Return a uniform floating point number on the interval 0<=x<1 */
375real
376gmx_rng_uniform_real(gmx_rng_t rng)
377{
378 if (sizeof(real) == sizeof(double))
379 {
380 return ((double)gmx_rng_uniform_uint32(rng))*(1.0/4294967296.0);
381 }
382 else
383 {
384 return ((float)gmx_rng_uniform_uint32(rng))*(1.0/4294967423.0);
385 }
386 /* divided by the smallest number that will generate a
387 * single precision real number on 0<=x<1.
388 * This needs to be slightly larger than MAX_UNIT since
389 * we are limited to an accuracy of 1e-7.
390 */
391}
392real
393
394gmx_rng_gaussian_table(gmx_rng_t rng)
395{
396 unsigned int i;
397
398 i = gmx_rng_uniform_uint32(rng);
399
400 /* The Gaussian table is a static constant in this file */
401 return gaussian_table[i >> (32 - GAUSS_TABLE14)];
402}
403
404void
405gmx_rng_cycle_2uniform(gmx_int64_t ctr1, gmx_int64_t ctr2,
406 gmx_int64_t key1, gmx_int64_t key2,
407 double* rnd)
408{
409 const gmx_int64_t mask_53bits = 0x1FFFFFFFFFFFFF;
410 const double two_power_min53 = 1.0/9007199254740992.0;
411
412 threefry2x64_ctr_t ctr = {{ctr1, ctr2}};
413 threefry2x64_key_t key = {{key1, key2}};
414 threefry2x64_ctr_t rand = threefry2x64(ctr, key)threefry2x64_R(threefry2x64_rounds, ctr, key);
415
416 rnd[0] = (rand.v[0] & mask_53bits)*two_power_min53;
417 rnd[1] = (rand.v[1] & mask_53bits)*two_power_min53;
418}
419
420void
421gmx_rng_cycle_3gaussian_table(gmx_int64_t ctr1, gmx_int64_t ctr2,
422 gmx_int64_t key1, gmx_int64_t key2,
423 real* rnd)
424{
425 threefry2x64_ctr_t ctr = {{ctr1, ctr2}};
426 threefry2x64_key_t key = {{key1, key2}};
427 threefry2x64_ctr_t rand = threefry2x64(ctr, key)threefry2x64_R(threefry2x64_rounds, ctr, key);
428
429 rnd[0] = gaussian_table[(rand.v[0] >> 48) & GAUSS_MASK((1 << 14) - 1)];
430 rnd[1] = gaussian_table[(rand.v[0] >> 32) & GAUSS_MASK((1 << 14) - 1)];
431 rnd[2] = gaussian_table[(rand.v[0] >> 16) & GAUSS_MASK((1 << 14) - 1)];
432}
433
434void
435gmx_rng_cycle_6gaussian_table(gmx_int64_t ctr1, gmx_int64_t ctr2,
436 gmx_int64_t key1, gmx_int64_t key2,
437 real* rnd)
438{
439 threefry2x64_ctr_t ctr = {{ctr1, ctr2}};
440 threefry2x64_key_t key = {{key1, key2}};
441 threefry2x64_ctr_t rand = threefry2x64(ctr, key)threefry2x64_R(threefry2x64_rounds, ctr, key);
442
443 rnd[0] = gaussian_table[(rand.v[0] >> 48) & GAUSS_MASK((1 << 14) - 1)];
444 rnd[1] = gaussian_table[(rand.v[0] >> 32) & GAUSS_MASK((1 << 14) - 1)];
445 rnd[2] = gaussian_table[(rand.v[0] >> 16) & GAUSS_MASK((1 << 14) - 1)];
446 rnd[3] = gaussian_table[(rand.v[1] >> 48) & GAUSS_MASK((1 << 14) - 1)];
447 rnd[4] = gaussian_table[(rand.v[1] >> 32) & GAUSS_MASK((1 << 14) - 1)];
448 rnd[5] = gaussian_table[(rand.v[1] >> 16) & GAUSS_MASK((1 << 14) - 1)];
449}