1/* GLIB - Library of useful routines for C programming
2 * Copyright (C) 1995-1997 Peter Mattis, Spencer Kimball and Josh MacDonald
3 *
4 * This library is free software; you can redistribute it and/or
5 * modify it under the terms of the GNU Lesser General Public
6 * License as published by the Free Software Foundation; either
7 * version 2.1 of the License, or (at your option) any later version.
8 *
9 * This library is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 * Lesser General Public License for more details.
13 *
14 * You should have received a copy of the GNU Lesser General Public
15 * License along with this library; if not, see <http://www.gnu.org/licenses/>.
16 */
17
18/* Originally developed and coded by Makoto Matsumoto and Takuji
19 * Nishimura. Please mail <matumoto@math.keio.ac.jp>, if you're using
20 * code from this file in your own programs or libraries.
21 * Further information on the Mersenne Twister can be found at
22 * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
23 * This code was adapted to glib by Sebastian Wilhelmi.
24 */
25
26/*
27 * Modified by the GLib Team and others 1997-2000. See the AUTHORS
28 * file for a list of people on the GLib Team. See the ChangeLog
29 * files for a list of changes. These files are distributed with
30 * GLib at ftp://ftp.gtk.org/pub/gtk/.
31 */
32
33/*
34 * MT safe
35 */
36
37#include "config.h"
38#define _CRT_RAND_S
39
40#include <math.h>
41#include <errno.h>
42#include <stdio.h>
43#include <string.h>
44#include <sys/types.h>
45#include "grand.h"
46
47#include "genviron.h"
48#include "gmain.h"
49#include "gmem.h"
50#include "gtestutils.h"
51#include "gthread.h"
52#include "gtimer.h"
53
54#ifdef G_OS_UNIX
55#include <unistd.h>
56#endif
57
58#ifdef G_OS_WIN32
59#include <stdlib.h>
60#include <process.h> /* For getpid() */
61#endif
62
63/**
64 * SECTION:random_numbers
65 * @title: Random Numbers
66 * @short_description: pseudo-random number generator
67 *
68 * The following functions allow you to use a portable, fast and good
69 * pseudo-random number generator (PRNG).
70 *
71 * Do not use this API for cryptographic purposes such as key
72 * generation, nonces, salts or one-time pads.
73 *
74 * This PRNG is suitable for non-cryptographic use such as in games
75 * (shuffling a card deck, generating levels), generating data for
76 * a test suite, etc. If you need random data for cryptographic
77 * purposes, it is recommended to use platform-specific APIs such
78 * as `/dev/random` on UNIX, or CryptGenRandom() on Windows.
79 *
80 * GRand uses the Mersenne Twister PRNG, which was originally
81 * developed by Makoto Matsumoto and Takuji Nishimura. Further
82 * information can be found at
83 * [this page](http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html).
84 *
85 * If you just need a random number, you simply call the g_random_*
86 * functions, which will create a globally used #GRand and use the
87 * according g_rand_* functions internally. Whenever you need a
88 * stream of reproducible random numbers, you better create a
89 * #GRand yourself and use the g_rand_* functions directly, which
90 * will also be slightly faster. Initializing a #GRand with a
91 * certain seed will produce exactly the same series of random
92 * numbers on all platforms. This can thus be used as a seed for
93 * e.g. games.
94 *
95 * The g_rand*_range functions will return high quality equally
96 * distributed random numbers, whereas for example the
97 * `(g_random_int()%max)` approach often
98 * doesn't yield equally distributed numbers.
99 *
100 * GLib changed the seeding algorithm for the pseudo-random number
101 * generator Mersenne Twister, as used by #GRand. This was necessary,
102 * because some seeds would yield very bad pseudo-random streams.
103 * Also the pseudo-random integers generated by g_rand*_int_range()
104 * will have a slightly better equal distribution with the new
105 * version of GLib.
106 *
107 * The original seeding and generation algorithms, as found in
108 * GLib 2.0.x, can be used instead of the new ones by setting the
109 * environment variable `G_RANDOM_VERSION` to the value of '2.0'.
110 * Use the GLib-2.0 algorithms only if you have sequences of numbers
111 * generated with Glib-2.0 that you need to reproduce exactly.
112 */
113
114/**
115 * GRand:
116 *
117 * The GRand struct is an opaque data structure. It should only be
118 * accessed through the g_rand_* functions.
119 **/
120
121G_LOCK_DEFINE_STATIC (global_random);
122
123/* Period parameters */
124#define N 624
125#define M 397
126#define MATRIX_A 0x9908b0df /* constant vector a */
127#define UPPER_MASK 0x80000000 /* most significant w-r bits */
128#define LOWER_MASK 0x7fffffff /* least significant r bits */
129
130/* Tempering parameters */
131#define TEMPERING_MASK_B 0x9d2c5680
132#define TEMPERING_MASK_C 0xefc60000
133#define TEMPERING_SHIFT_U(y) (y >> 11)
134#define TEMPERING_SHIFT_S(y) (y << 7)
135#define TEMPERING_SHIFT_T(y) (y << 15)
136#define TEMPERING_SHIFT_L(y) (y >> 18)
137
138static guint
139get_random_version (void)
140{
141 static gsize initialized = FALSE;
142 static guint random_version;
143
144 if (g_once_init_enter (&initialized))
145 {
146 const gchar *version_string = g_getenv (variable: "G_RANDOM_VERSION");
147 if (!version_string || version_string[0] == '\000' ||
148 strcmp (s1: version_string, s2: "2.2") == 0)
149 random_version = 22;
150 else if (strcmp (s1: version_string, s2: "2.0") == 0)
151 random_version = 20;
152 else
153 {
154 g_warning ("Unknown G_RANDOM_VERSION \"%s\". Using version 2.2.",
155 version_string);
156 random_version = 22;
157 }
158 g_once_init_leave (&initialized, TRUE);
159 }
160
161 return random_version;
162}
163
164struct _GRand
165{
166 guint32 mt[N]; /* the array for the state vector */
167 guint mti;
168};
169
170/**
171 * g_rand_new_with_seed:
172 * @seed: a value to initialize the random number generator
173 *
174 * Creates a new random number generator initialized with @seed.
175 *
176 * Returns: the new #GRand
177 **/
178GRand*
179g_rand_new_with_seed (guint32 seed)
180{
181 GRand *rand = g_new0 (GRand, 1);
182 g_rand_set_seed (rand_: rand, seed);
183 return rand;
184}
185
186/**
187 * g_rand_new_with_seed_array:
188 * @seed: an array of seeds to initialize the random number generator
189 * @seed_length: an array of seeds to initialize the random number
190 * generator
191 *
192 * Creates a new random number generator initialized with @seed.
193 *
194 * Returns: the new #GRand
195 *
196 * Since: 2.4
197 */
198GRand*
199g_rand_new_with_seed_array (const guint32 *seed,
200 guint seed_length)
201{
202 GRand *rand = g_new0 (GRand, 1);
203 g_rand_set_seed_array (rand_: rand, seed, seed_length);
204 return rand;
205}
206
207/**
208 * g_rand_new:
209 *
210 * Creates a new random number generator initialized with a seed taken
211 * either from `/dev/urandom` (if existing) or from the current time
212 * (as a fallback).
213 *
214 * On Windows, the seed is taken from rand_s().
215 *
216 * Returns: the new #GRand
217 */
218GRand*
219g_rand_new (void)
220{
221 guint32 seed[4];
222#ifdef G_OS_UNIX
223 static gboolean dev_urandom_exists = TRUE;
224
225 if (dev_urandom_exists)
226 {
227 FILE* dev_urandom;
228
229 do
230 {
231 dev_urandom = fopen(filename: "/dev/urandom", modes: "rb");
232 }
233 while G_UNLIKELY (dev_urandom == NULL && errno == EINTR);
234
235 if (dev_urandom)
236 {
237 int r;
238
239 setvbuf (stream: dev_urandom, NULL, _IONBF, n: 0);
240 do
241 {
242 errno = 0;
243 r = fread (ptr: seed, size: sizeof (seed), n: 1, stream: dev_urandom);
244 }
245 while G_UNLIKELY (errno == EINTR);
246
247 if (r != 1)
248 dev_urandom_exists = FALSE;
249
250 fclose (stream: dev_urandom);
251 }
252 else
253 dev_urandom_exists = FALSE;
254 }
255
256 if (!dev_urandom_exists)
257 {
258 gint64 now_us = g_get_real_time ();
259 seed[0] = now_us / G_USEC_PER_SEC;
260 seed[1] = now_us % G_USEC_PER_SEC;
261 seed[2] = getpid ();
262 seed[3] = getppid ();
263 }
264#else /* G_OS_WIN32 */
265 /* rand_s() is only available since Visual Studio 2005 and
266 * MinGW-w64 has a wrapper that will emulate rand_s() if it's not in msvcrt
267 */
268#if (defined(_MSC_VER) && _MSC_VER >= 1400) || defined(__MINGW64_VERSION_MAJOR)
269 gint i;
270
271 for (i = 0; i < G_N_ELEMENTS (seed); i++)
272 rand_s (&seed[i]);
273#else
274#warning Using insecure seed for random number generation because of missing rand_s() in Windows XP
275 GTimeVal now;
276
277 g_get_current_time (&now);
278 seed[0] = now.tv_sec;
279 seed[1] = now.tv_usec;
280 seed[2] = getpid ();
281 seed[3] = 0;
282#endif
283
284#endif
285
286 return g_rand_new_with_seed_array (seed, seed_length: 4);
287}
288
289/**
290 * g_rand_free:
291 * @rand_: a #GRand
292 *
293 * Frees the memory allocated for the #GRand.
294 */
295void
296g_rand_free (GRand *rand)
297{
298 g_return_if_fail (rand != NULL);
299
300 g_free (mem: rand);
301}
302
303/**
304 * g_rand_copy:
305 * @rand_: a #GRand
306 *
307 * Copies a #GRand into a new one with the same exact state as before.
308 * This way you can take a snapshot of the random number generator for
309 * replaying later.
310 *
311 * Returns: the new #GRand
312 *
313 * Since: 2.4
314 */
315GRand*
316g_rand_copy (GRand *rand)
317{
318 GRand* new_rand;
319
320 g_return_val_if_fail (rand != NULL, NULL);
321
322 new_rand = g_new0 (GRand, 1);
323 memcpy (dest: new_rand, src: rand, n: sizeof (GRand));
324
325 return new_rand;
326}
327
328/**
329 * g_rand_set_seed:
330 * @rand_: a #GRand
331 * @seed: a value to reinitialize the random number generator
332 *
333 * Sets the seed for the random number generator #GRand to @seed.
334 */
335void
336g_rand_set_seed (GRand *rand,
337 guint32 seed)
338{
339 g_return_if_fail (rand != NULL);
340
341 switch (get_random_version ())
342 {
343 case 20:
344 /* setting initial seeds to mt[N] using */
345 /* the generator Line 25 of Table 1 in */
346 /* [KNUTH 1981, The Art of Computer Programming */
347 /* Vol. 2 (2nd Ed.), pp102] */
348
349 if (seed == 0) /* This would make the PRNG produce only zeros */
350 seed = 0x6b842128; /* Just set it to another number */
351
352 rand->mt[0]= seed;
353 for (rand->mti=1; rand->mti<N; rand->mti++)
354 rand->mt[rand->mti] = (69069 * rand->mt[rand->mti-1]);
355
356 break;
357 case 22:
358 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
359 /* In the previous version (see above), MSBs of the */
360 /* seed affect only MSBs of the array mt[]. */
361
362 rand->mt[0]= seed;
363 for (rand->mti=1; rand->mti<N; rand->mti++)
364 rand->mt[rand->mti] = 1812433253UL *
365 (rand->mt[rand->mti-1] ^ (rand->mt[rand->mti-1] >> 30)) + rand->mti;
366 break;
367 default:
368 g_assert_not_reached ();
369 }
370}
371
372/**
373 * g_rand_set_seed_array:
374 * @rand_: a #GRand
375 * @seed: array to initialize with
376 * @seed_length: length of array
377 *
378 * Initializes the random number generator by an array of longs.
379 * Array can be of arbitrary size, though only the first 624 values
380 * are taken. This function is useful if you have many low entropy
381 * seeds, or if you require more then 32 bits of actual entropy for
382 * your application.
383 *
384 * Since: 2.4
385 */
386void
387g_rand_set_seed_array (GRand *rand,
388 const guint32 *seed,
389 guint seed_length)
390{
391 guint i, j, k;
392
393 g_return_if_fail (rand != NULL);
394 g_return_if_fail (seed_length >= 1);
395
396 g_rand_set_seed (rand, seed: 19650218UL);
397
398 i=1; j=0;
399 k = (N>seed_length ? N : seed_length);
400 for (; k; k--)
401 {
402 rand->mt[i] = (rand->mt[i] ^
403 ((rand->mt[i-1] ^ (rand->mt[i-1] >> 30)) * 1664525UL))
404 + seed[j] + j; /* non linear */
405 rand->mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
406 i++; j++;
407 if (i>=N)
408 {
409 rand->mt[0] = rand->mt[N-1];
410 i=1;
411 }
412 if (j>=seed_length)
413 j=0;
414 }
415 for (k=N-1; k; k--)
416 {
417 rand->mt[i] = (rand->mt[i] ^
418 ((rand->mt[i-1] ^ (rand->mt[i-1] >> 30)) * 1566083941UL))
419 - i; /* non linear */
420 rand->mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
421 i++;
422 if (i>=N)
423 {
424 rand->mt[0] = rand->mt[N-1];
425 i=1;
426 }
427 }
428
429 rand->mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
430}
431
432/**
433 * g_rand_boolean:
434 * @rand_: a #GRand
435 *
436 * Returns a random #gboolean from @rand_.
437 * This corresponds to an unbiased coin toss.
438 *
439 * Returns: a random #gboolean
440 */
441/**
442 * g_rand_int:
443 * @rand_: a #GRand
444 *
445 * Returns the next random #guint32 from @rand_ equally distributed over
446 * the range [0..2^32-1].
447 *
448 * Returns: a random number
449 */
450guint32
451g_rand_int (GRand *rand)
452{
453 guint32 y;
454 static const guint32 mag01[2]={0x0, MATRIX_A};
455 /* mag01[x] = x * MATRIX_A for x=0,1 */
456
457 g_return_val_if_fail (rand != NULL, 0);
458
459 if (rand->mti >= N) { /* generate N words at one time */
460 int kk;
461
462 for (kk = 0; kk < N - M; kk++) {
463 y = (rand->mt[kk]&UPPER_MASK)|(rand->mt[kk+1]&LOWER_MASK);
464 rand->mt[kk] = rand->mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1];
465 }
466 for (; kk < N - 1; kk++) {
467 y = (rand->mt[kk]&UPPER_MASK)|(rand->mt[kk+1]&LOWER_MASK);
468 rand->mt[kk] = rand->mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1];
469 }
470 y = (rand->mt[N-1]&UPPER_MASK)|(rand->mt[0]&LOWER_MASK);
471 rand->mt[N-1] = rand->mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1];
472
473 rand->mti = 0;
474 }
475
476 y = rand->mt[rand->mti++];
477 y ^= TEMPERING_SHIFT_U(y);
478 y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
479 y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
480 y ^= TEMPERING_SHIFT_L(y);
481
482 return y;
483}
484
485/* transform [0..2^32] -> [0..1] */
486#define G_RAND_DOUBLE_TRANSFORM 2.3283064365386962890625e-10
487
488/**
489 * g_rand_int_range:
490 * @rand_: a #GRand
491 * @begin: lower closed bound of the interval
492 * @end: upper open bound of the interval
493 *
494 * Returns the next random #gint32 from @rand_ equally distributed over
495 * the range [@begin..@end-1].
496 *
497 * Returns: a random number
498 */
499gint32
500g_rand_int_range (GRand *rand,
501 gint32 begin,
502 gint32 end)
503{
504 guint32 dist = end - begin;
505 guint32 random = 0;
506
507 g_return_val_if_fail (rand != NULL, begin);
508 g_return_val_if_fail (end > begin, begin);
509
510 switch (get_random_version ())
511 {
512 case 20:
513 if (dist <= 0x10000L) /* 2^16 */
514 {
515 /* This method, which only calls g_rand_int once is only good
516 * for (end - begin) <= 2^16, because we only have 32 bits set
517 * from the one call to g_rand_int ().
518 *
519 * We are using (trans + trans * trans), because g_rand_int only
520 * covers [0..2^32-1] and thus g_rand_int * trans only covers
521 * [0..1-2^-32], but the biggest double < 1 is 1-2^-52.
522 */
523
524 gdouble double_rand = g_rand_int (rand) *
525 (G_RAND_DOUBLE_TRANSFORM +
526 G_RAND_DOUBLE_TRANSFORM * G_RAND_DOUBLE_TRANSFORM);
527
528 random = (gint32) (double_rand * dist);
529 }
530 else
531 {
532 /* Now we use g_rand_double_range (), which will set 52 bits
533 * for us, so that it is safe to round and still get a decent
534 * distribution
535 */
536 random = (gint32) g_rand_double_range (rand_: rand, begin: 0, end: dist);
537 }
538 break;
539 case 22:
540 if (dist == 0)
541 random = 0;
542 else
543 {
544 /* maxvalue is set to the predecessor of the greatest
545 * multiple of dist less or equal 2^32.
546 */
547 guint32 maxvalue;
548 if (dist <= 0x80000000u) /* 2^31 */
549 {
550 /* maxvalue = 2^32 - 1 - (2^32 % dist) */
551 guint32 leftover = (0x80000000u % dist) * 2;
552 if (leftover >= dist) leftover -= dist;
553 maxvalue = 0xffffffffu - leftover;
554 }
555 else
556 maxvalue = dist - 1;
557
558 do
559 random = g_rand_int (rand);
560 while (random > maxvalue);
561
562 random %= dist;
563 }
564 break;
565 default:
566 g_assert_not_reached ();
567 }
568
569 return begin + random;
570}
571
572/**
573 * g_rand_double:
574 * @rand_: a #GRand
575 *
576 * Returns the next random #gdouble from @rand_ equally distributed over
577 * the range [0..1).
578 *
579 * Returns: a random number
580 */
581gdouble
582g_rand_double (GRand *rand)
583{
584 /* We set all 52 bits after the point for this, not only the first
585 32. That's why we need two calls to g_rand_int */
586 gdouble retval = g_rand_int (rand) * G_RAND_DOUBLE_TRANSFORM;
587 retval = (retval + g_rand_int (rand)) * G_RAND_DOUBLE_TRANSFORM;
588
589 /* The following might happen due to very bad rounding luck, but
590 * actually this should be more than rare, we just try again then */
591 if (retval >= 1.0)
592 return g_rand_double (rand);
593
594 return retval;
595}
596
597/**
598 * g_rand_double_range:
599 * @rand_: a #GRand
600 * @begin: lower closed bound of the interval
601 * @end: upper open bound of the interval
602 *
603 * Returns the next random #gdouble from @rand_ equally distributed over
604 * the range [@begin..@end).
605 *
606 * Returns: a random number
607 */
608gdouble
609g_rand_double_range (GRand *rand,
610 gdouble begin,
611 gdouble end)
612{
613 gdouble r;
614
615 r = g_rand_double (rand);
616
617 return r * end - (r - 1) * begin;
618}
619
620static GRand *
621get_global_random (void)
622{
623 static GRand *global_random;
624
625 /* called while locked */
626 if (!global_random)
627 global_random = g_rand_new ();
628
629 return global_random;
630}
631
632/**
633 * g_random_boolean:
634 *
635 * Returns a random #gboolean.
636 * This corresponds to an unbiased coin toss.
637 *
638 * Returns: a random #gboolean
639 */
640/**
641 * g_random_int:
642 *
643 * Return a random #guint32 equally distributed over the range
644 * [0..2^32-1].
645 *
646 * Returns: a random number
647 */
648guint32
649g_random_int (void)
650{
651 guint32 result;
652 G_LOCK (global_random);
653 result = g_rand_int (rand: get_global_random ());
654 G_UNLOCK (global_random);
655 return result;
656}
657
658/**
659 * g_random_int_range:
660 * @begin: lower closed bound of the interval
661 * @end: upper open bound of the interval
662 *
663 * Returns a random #gint32 equally distributed over the range
664 * [@begin..@end-1].
665 *
666 * Returns: a random number
667 */
668gint32
669g_random_int_range (gint32 begin,
670 gint32 end)
671{
672 gint32 result;
673 G_LOCK (global_random);
674 result = g_rand_int_range (rand: get_global_random (), begin, end);
675 G_UNLOCK (global_random);
676 return result;
677}
678
679/**
680 * g_random_double:
681 *
682 * Returns a random #gdouble equally distributed over the range [0..1).
683 *
684 * Returns: a random number
685 */
686gdouble
687g_random_double (void)
688{
689 double result;
690 G_LOCK (global_random);
691 result = g_rand_double (rand: get_global_random ());
692 G_UNLOCK (global_random);
693 return result;
694}
695
696/**
697 * g_random_double_range:
698 * @begin: lower closed bound of the interval
699 * @end: upper open bound of the interval
700 *
701 * Returns a random #gdouble equally distributed over the range
702 * [@begin..@end).
703 *
704 * Returns: a random number
705 */
706gdouble
707g_random_double_range (gdouble begin,
708 gdouble end)
709{
710 double result;
711 G_LOCK (global_random);
712 result = g_rand_double_range (rand: get_global_random (), begin, end);
713 G_UNLOCK (global_random);
714 return result;
715}
716
717/**
718 * g_random_set_seed:
719 * @seed: a value to reinitialize the global random number generator
720 *
721 * Sets the seed for the global random number generator, which is used
722 * by the g_random_* functions, to @seed.
723 */
724void
725g_random_set_seed (guint32 seed)
726{
727 G_LOCK (global_random);
728 g_rand_set_seed (rand: get_global_random (), seed);
729 G_UNLOCK (global_random);
730}
731

source code of gtk/subprojects/glib/glib/grand.c