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