glib/grand.c

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