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 "glib.h"
  51 #include "gthreadprivate.h"
  52 #include "galias.h"
  53
  54 #ifdef G_OS_WIN32
  55 #include <process.h>            /* For getpid() */
  56 #endif
  57
  58 G_LOCK_DEFINE_STATIC (global_random);
  59 static GRand* global_random = NULL;
  60
  61 /* Period parameters */
  62 #define N 624
  63 #define M 397
  64 #define MATRIX_A 0x9908b0df   /* constant vector a */
  65 #define UPPER_MASK 0x80000000 /* most significant w-r bits */
  66 #define LOWER_MASK 0x7fffffff /* least significant r bits */
  67
  68 /* Tempering parameters */
  69 #define TEMPERING_MASK_B 0x9d2c5680
  70 #define TEMPERING_MASK_C 0xefc60000
  71 #define TEMPERING_SHIFT_U(y)  (y >> 11)
  72 #define TEMPERING_SHIFT_S(y)  (y << 7)
  73 #define TEMPERING_SHIFT_T(y)  (y << 15)
  74 #define TEMPERING_SHIFT_L(y)  (y >> 18)
  75
  76 static guint
  77 get_random_version (void)
  78 {
  79   static gboolean initialized = FALSE;
  80   static guint random_version;
  81
  82   if (!initialized)
  83     {
  84       const gchar *version_string = g_getenv ("G_RANDOM_VERSION");
  85       if (!version_string || version_string[0] == '\000' ||
  86           strcmp (version_string, "2.2") == 0)
  87         random_version = 22;
  88       else if (strcmp (version_string, "2.0") == 0)
  89         random_version = 20;
  90       else
  91         {
  92           g_warning ("Unknown G_RANDOM_VERSION \"%s\". Using version 2.2.",
  93                      version_string);
  94           random_version = 22;
  95         }
  96       initialized = TRUE;
  97     }
  98
  99   return random_version;
 100 }
 101
 102 /* This is called from g_thread_init(). It's used to
 103  * initialize some static data in a threadsafe way.
 104  */
 105 void
 106 _g_rand_thread_init (void)
 107 {
 108   (void)get_random_version ();
 109 }
 110
 111 struct _GRand
 112 {
 113   guint32 mt[N]; /* the array for the state vector  */
 114   guint mti;
 115 };
 116
 117 /**
 118  * g_rand_new_with_seed:
 119  * @seed: a value to initialize the random number generator.
 120  *
 121  * Creates a new random number generator initialized with @seed.
 122  *
 123  * Return value: the new #GRand.
 124  **/
 125 GRand*
 126 g_rand_new_with_seed (guint32 seed)
 127 {
 128   GRand *rand = g_new0 (GRand, 1);
 129   g_rand_set_seed (rand, seed);
 130   return rand;
 131 }
 132
 133 /**
 134  * g_rand_new_with_seed_array:
 135  * @seed: an array of seeds to initialize the random number generator.
 136  * @seed_length: an array of seeds to initialize the random number generator.
 137  *
 138  * Creates a new random number generator initialized with @seed.
 139  *
 140  * Return value: the new #GRand.
 141  *
 142  * Since: 2.4
 143  **/
 144 GRand*
 145 g_rand_new_with_seed_array (const guint32 *seed, guint seed_length)
 146 {
 147   GRand *rand = g_new0 (GRand, 1);
 148   g_rand_set_seed_array (rand, seed, seed_length);
 149   return rand;
 150 }
 151
 152 /**
 153  * g_rand_new:
 154  *
 155  * Creates a new random number generator initialized with a seed taken
 156  * either from <filename>/dev/urandom</filename> (if existing) or from
 157  * the current time (as a fallback).
 158  *
 159  * Return value: the new #GRand.
 160  **/
 161 GRand*
 162 g_rand_new (void)
 163 {
 164   guint32 seed[4];
 165   GTimeVal now;
 166 #ifdef G_OS_UNIX
 167   static gboolean dev_urandom_exists = TRUE;
 168
 169   if (dev_urandom_exists)
 170     {
 171       FILE* dev_urandom;
 172
 173       do
 174         {
 175           errno = 0;
 176           dev_urandom = fopen("/dev/urandom", "rb");
 177         }
 178       while G_UNLIKELY (errno == EINTR);
 179
 180       if (dev_urandom)
 181         {
 182           int r;
 183
 184           setvbuf (dev_urandom, NULL, _IONBF, 0);
 185           do
 186             {
 187               errno = 0;
 188               r = fread (seed, sizeof (seed), 1, dev_urandom);
 189             }
 190           while G_UNLIKELY (errno == EINTR);
 191
 192           if (r != 1)
 193             dev_urandom_exists = FALSE;
 194
 195           fclose (dev_urandom);
 196         }
 197       else
 198         dev_urandom_exists = FALSE;
 199     }
 200 #else
 201   static gboolean dev_urandom_exists = FALSE;
 202 #endif
 203
 204   if (!dev_urandom_exists)
 205     {
 206       g_get_current_time (&now);
 207       seed[0] = now.tv_sec;
 208       seed[1] = now.tv_usec;
 209       seed[2] = getpid ();
 210 #ifdef G_OS_UNIX
 211       seed[3] = getppid ();
 212 #else
 213       seed[3] = 0;
 214 #endif
 215     }
 216
 217   return g_rand_new_with_seed_array (seed, 4);
 218 }
 219
 220 /**
 221  * g_rand_free:
 222  * @rand_: a #GRand.
 223  *
 224  * Frees the memory allocated for the #GRand.
 225  **/
 226 void
 227 g_rand_free (GRand* rand)
 228 {
 229   g_return_if_fail (rand != NULL);
 230
 231   g_free (rand);
 232 }
 233
 234 /**
 235  * g_rand_copy:
 236  * @rand_: a #GRand.
 237  *
 238  * Copies a #GRand into a new one with the same exact state as before.
 239  * This way you can take a snapshot of the random number generator for
 240  * replaying later.
 241  *
 242  * Return value: the new #GRand.
 243  *
 244  * Since: 2.4
 245  **/
 246 GRand *
 247 g_rand_copy (GRand* rand)
 248 {
 249   GRand* new_rand;
 250
 251   g_return_val_if_fail (rand != NULL, NULL);
 252
 253   new_rand = g_new0 (GRand, 1);
 254   memcpy (new_rand, rand, sizeof (GRand));
 255
 256   return new_rand;
 257 }
 258
 259 /**
 260  * g_rand_set_seed:
 261  * @rand_: a #GRand.
 262  * @seed: a value to reinitialize the random number generator.
 263  *
 264  * Sets the seed for the random number generator #GRand to @seed.
 265  **/
 266 void
 267 g_rand_set_seed (GRand* rand, guint32 seed)
 268 {
 269   g_return_if_fail (rand != NULL);
 270
 271   switch (get_random_version ())
 272     {
 273     case 20:
 274       /* setting initial seeds to mt[N] using         */
 275       /* the generator Line 25 of Table 1 in          */
 276       /* [KNUTH 1981, The Art of Computer Programming */
 277       /*    Vol. 2 (2nd Ed.), pp102]                  */
 278
 279       if (seed == 0) /* This would make the PRNG procude only zeros */
 280         seed = 0x6b842128; /* Just set it to another number */
 281
 282       rand->mt[0]= seed;
 283       for (rand->mti=1; rand->mti<N; rand->mti++)
 284         rand->mt[rand->mti] = (69069 * rand->mt[rand->mti-1]);
 285
 286       break;
 287     case 22:
 288       /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
 289       /* In the previous version (see above), MSBs of the    */
 290       /* seed affect only MSBs of the array mt[].            */
 291
 292       rand->mt[0]= seed;
 293       for (rand->mti=1; rand->mti<N; rand->mti++)
 294         rand->mt[rand->mti] = 1812433253UL *
 295           (rand->mt[rand->mti-1] ^ (rand->mt[rand->mti-1] >> 30)) + rand->mti;
 296       break;
 297     default:
 298       g_assert_not_reached ();
 299     }
 300 }
 301
 302 /**
 303  * g_rand_set_seed_array:
 304  * @rand_: a #GRand.
 305  * @seed: array to initialize with
 306  * @seed_length: length of array
 307  *
 308  * Initializes the random number generator by an array of
 309  * longs.  Array can be of arbitrary size, though only the
 310  * first 624 values are taken.  This function is useful
 311  * if you have many low entropy seeds, or if you require more then
 312  * 32bits of actual entropy for your application.
 313  *
 314  * Since: 2.4
 315  **/
 316 void
 317 g_rand_set_seed_array (GRand* rand, const guint32 *seed, guint seed_length)
 318 {
 319   int i, j, k;
 320
 321   g_return_if_fail (rand != NULL);
 322   g_return_if_fail (seed_length >= 1);
 323
 324   g_rand_set_seed (rand, 19650218UL);
 325
 326   i=1; j=0;
 327   k = (N>seed_length ? N : seed_length);
 328   for (; k; k--)
 329     {
 330       rand->mt[i] = (rand->mt[i] ^
 331                      ((rand->mt[i-1] ^ (rand->mt[i-1] >> 30)) * 1664525UL))
 332               + seed[j] + j; /* non linear */
 333       rand->mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
 334       i++; j++;
 335       if (i>=N)
 336         {
 337           rand->mt[0] = rand->mt[N-1];
 338           i=1;
 339         }
 340       if (j>=seed_length)
 341         j=0;
 342     }
 343   for (k=N-1; k; k--)
 344     {
 345       rand->mt[i] = (rand->mt[i] ^
 346                      ((rand->mt[i-1] ^ (rand->mt[i-1] >> 30)) * 1566083941UL))
 347               - i; /* non linear */
 348       rand->mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
 349       i++;
 350       if (i>=N)
 351         {
 352           rand->mt[0] = rand->mt[N-1];
 353           i=1;
 354         }
 355     }
 356
 357   rand->mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
 358 }
 359
 360 /**
 361  * g_rand_int:
 362  * @rand_: a #GRand.
 363  *
 364  * Returns the next random #guint32 from @rand_ equally distributed over
 365  * the range [0..2^32-1].
 366  *
 367  * Return value: A random number.
 368  **/
 369 guint32
 370 g_rand_int (GRand* rand)
 371 {
 372   guint32 y;
 373   static const guint32 mag01[2]={0x0, MATRIX_A};
 374   /* mag01[x] = x * MATRIX_A  for x=0,1 */
 375
 376   g_return_val_if_fail (rand != NULL, 0);
 377
 378   if (rand->mti >= N) { /* generate N words at one time */
 379     int kk;
 380
 381     for (kk=0;kk<N-M;kk++) {
 382       y = (rand->mt[kk]&UPPER_MASK)|(rand->mt[kk+1]&LOWER_MASK);
 383       rand->mt[kk] = rand->mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1];
 384     }
 385     for (;kk<N-1;kk++) {
 386       y = (rand->mt[kk]&UPPER_MASK)|(rand->mt[kk+1]&LOWER_MASK);
 387       rand->mt[kk] = rand->mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1];
 388     }
 389     y = (rand->mt[N-1]&UPPER_MASK)|(rand->mt[0]&LOWER_MASK);
 390     rand->mt[N-1] = rand->mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1];
 391
 392     rand->mti = 0;
 393   }
 394
 395   y = rand->mt[rand->mti++];
 396   y ^= TEMPERING_SHIFT_U(y);
 397   y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
 398   y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
 399   y ^= TEMPERING_SHIFT_L(y);
 400
 401   return y;
 402 }
 403
 404 /* transform [0..2^32] -> [0..1] */
 405 #define G_RAND_DOUBLE_TRANSFORM 2.3283064365386962890625e-10
 406
 407 /**
 408  * g_rand_int_range:
 409  * @rand_: a #GRand.
 410  * @begin: lower closed bound of the interval.
 411  * @end: upper open bound of the interval.
 412  *
 413  * Returns the next random #gint32 from @rand_ equally distributed over
 414  * the range [@begin..@end-1].
 415  *
 416  * Return value: A random number.
 417  **/
 418 gint32
 419 g_rand_int_range (GRand* rand, gint32 begin, gint32 end)
 420 {
 421   guint32 dist = end - begin;
 422   guint32 random;
 423
 424   g_return_val_if_fail (rand != NULL, begin);
 425   g_return_val_if_fail (end > begin, begin);
 426
 427   switch (get_random_version ())
 428     {
 429     case 20:
 430       if (dist <= 0x10000L) /* 2^16 */
 431         {
 432           /* This method, which only calls g_rand_int once is only good
 433            * for (end - begin) <= 2^16, because we only have 32 bits set
 434            * from the one call to g_rand_int (). */
 435
 436           /* we are using (trans + trans * trans), because g_rand_int only
 437            * covers [0..2^32-1] and thus g_rand_int * trans only covers
 438            * [0..1-2^-32], but the biggest double < 1 is 1-2^-52.
 439            */
 440
 441           gdouble double_rand = g_rand_int (rand) *
 442             (G_RAND_DOUBLE_TRANSFORM +
 443              G_RAND_DOUBLE_TRANSFORM * G_RAND_DOUBLE_TRANSFORM);
 444
 445           random = (gint32) (double_rand * dist);
 446         }
 447       else
 448         {
 449           /* Now we use g_rand_double_range (), which will set 52 bits for
 450              us, so that it is safe to round and still get a decent
 451              distribution */
 452           random = (gint32) g_rand_double_range (rand, 0, dist);
 453         }
 454       break;
 455     case 22:
 456       if (dist == 0)
 457         random = 0;
 458       else
 459         {
 460           /* maxvalue is set to the predecessor of the greatest
 461            * multiple of dist less or equal 2^32. */
 462           guint32 maxvalue;
 463           if (dist <= 0x80000000u) /* 2^31 */
 464             {
 465               /* maxvalue = 2^32 - 1 - (2^32 % dist) */
 466               guint32 leftover = (0x80000000u % dist) * 2;
 467               if (leftover >= dist) leftover -= dist;
 468               maxvalue = 0xffffffffu - leftover;
 469             }
 470           else
 471             maxvalue = dist - 1;
 472
 473           do
 474             random = g_rand_int (rand);
 475           while (random > maxvalue);
 476
 477           random %= dist;
 478         }
 479       break;
 480     default:
 481       random = 0;               /* Quiet GCC */
 482       g_assert_not_reached ();
 483     }
 484
 485   return begin + random;
 486 }
 487
 488 /**
 489  * g_rand_double:
 490  * @rand_: a #GRand.
 491  *
 492  * Returns the next random #gdouble from @rand_ equally distributed over
 493  * the range [0..1).
 494  *
 495  * Return value: A random number.
 496  **/
 497 gdouble
 498 g_rand_double (GRand* rand)
 499 {
 500   /* We set all 52 bits after the point for this, not only the first
 501      32. Thats why we need two calls to g_rand_int */
 502   gdouble retval = g_rand_int (rand) * G_RAND_DOUBLE_TRANSFORM;
 503   retval = (retval + g_rand_int (rand)) * G_RAND_DOUBLE_TRANSFORM;
 504
 505   /* The following might happen due to very bad rounding luck, but
 506    * actually this should be more than rare, we just try again then */
 507   if (retval >= 1.0)
 508     return g_rand_double (rand);
 509
 510   return retval;
 511 }
 512
 513 /**
 514  * g_rand_double_range:
 515  * @rand_: a #GRand.
 516  * @begin: lower closed bound of the interval.
 517  * @end: upper open bound of the interval.
 518  *
 519  * Returns the next random #gdouble from @rand_ equally distributed over
 520  * the range [@begin..@end).
 521  *
 522  * Return value: A random number.
 523  **/
 524 gdouble
 525 g_rand_double_range (GRand* rand, gdouble begin, gdouble end)
 526 {
 527   return g_rand_double (rand) * (end - begin) + begin;
 528 }
 529
 530 /**
 531  * g_random_int:
 532  *
 533  * Return a random #guint32 equally distributed over the range
 534  * [0..2^32-1].
 535  *
 536  * Return value: A random number.
 537  **/
 538 guint32
 539 g_random_int (void)
 540 {
 541   guint32 result;
 542   G_LOCK (global_random);
 543   if (!global_random)
 544     global_random = g_rand_new ();
 545
 546   result = g_rand_int (global_random);
 547   G_UNLOCK (global_random);
 548   return result;
 549 }
 550
 551 /**
 552  * g_random_int_range:
 553  * @begin: lower closed bound of the interval.
 554  * @end: upper open bound of the interval.
 555  *
 556  * Returns a random #gint32 equally distributed over the range
 557  * [@begin..@end-1].
 558  *
 559  * Return value: A random number.
 560  **/
 561 gint32
 562 g_random_int_range (gint32 begin, gint32 end)
 563 {
 564   gint32 result;
 565   G_LOCK (global_random);
 566   if (!global_random)
 567     global_random = g_rand_new ();
 568
 569   result = g_rand_int_range (global_random, begin, end);
 570   G_UNLOCK (global_random);
 571   return result;
 572 }
 573
 574 /**
 575  * g_random_double:
 576  *
 577  * Returns a random #gdouble equally distributed over the range [0..1).
 578  *
 579  * Return value: A random number.
 580  **/
 581 gdouble
 582 g_random_double (void)
 583 {
 584   double result;
 585   G_LOCK (global_random);
 586   if (!global_random)
 587     global_random = g_rand_new ();
 588
 589   result = g_rand_double (global_random);
 590   G_UNLOCK (global_random);
 591   return result;
 592 }
 593
 594 /**
 595  * g_random_double_range:
 596  * @begin: lower closed bound of the interval.
 597  * @end: upper open bound of the interval.
 598  *
 599  * Returns a random #gdouble equally distributed over the range [@begin..@end).
 600  *
 601  * Return value: A random number.
 602  **/
 603 gdouble
 604 g_random_double_range (gdouble begin, gdouble end)
 605 {
 606   double result;
 607   G_LOCK (global_random);
 608   if (!global_random)
 609     global_random = g_rand_new ();
 610
 611   result = g_rand_double_range (global_random, begin, end);
 612   G_UNLOCK (global_random);
 613   return result;
 614 }
 615
 616 /**
 617  * g_random_set_seed:
 618  * @seed: a value to reinitialize the global random number generator.
 619  *
 620  * Sets the seed for the global random number generator, which is used
 621  * by the <function>g_random_*</function> functions, to @seed.
 622  **/
 623 void
 624 g_random_set_seed (guint32 seed)
 625 {
 626   G_LOCK (global_random);
 627   if (!global_random)
 628     global_random = g_rand_new_with_seed (seed);
 629   else
 630     g_rand_set_seed (global_random, seed);
 631   G_UNLOCK (global_random);
 632 }
 633
 634
 635 #define __G_RAND_C__
 636 #include "galiasdef.c"