From 95aff22dff87bd2c9db875a938e1e116ebdaf33c Mon Sep 17 00:00:00 2001 From: Sebastian Wilhelmi Date: Fri, 9 Apr 1999 14:40:58 +0000 Subject: [PATCH] New files to implement the Mersenne Twister Pseudo Random Number 1999-04-09 Sebastian Wilhelmi * grand.c, tests/rand-test.c: New files to implement the Mersenne Twister Pseudo Random Number Generator. * glib.h, AUTHORS, Makefile.am, tests/Makefile.am: Changed accordingly. --- AUTHORS | 5 + ChangeLog | 8 ++ ChangeLog.pre-2-0 | 8 ++ ChangeLog.pre-2-10 | 8 ++ ChangeLog.pre-2-12 | 8 ++ ChangeLog.pre-2-2 | 8 ++ ChangeLog.pre-2-4 | 8 ++ ChangeLog.pre-2-6 | 8 ++ ChangeLog.pre-2-8 | 8 ++ Makefile.am | 1 + glib.h | 43 +++++++ glib/Makefile.am | 1 + glib/glib.h | 43 +++++++ glib/grand.c | 334 ++++++++++++++++++++++++++++++++++++++++++++++++++++ grand.c | 334 ++++++++++++++++++++++++++++++++++++++++++++++++++++ gthread/Makefile.am | 2 +- tests/Makefile.am | 2 + tests/rand-test.c | 43 +++++++ 18 files changed, 871 insertions(+), 1 deletion(-) create mode 100644 glib/grand.c create mode 100644 grand.c create mode 100644 tests/rand-test.c diff --git a/AUTHORS b/AUTHORS index e898236..49c3037 100644 --- a/AUTHORS +++ b/AUTHORS @@ -23,3 +23,8 @@ Sebastian Wilhelmi There are also many others who have contributed patches and fixes; we thank them, for helping us in advancing GLIB. + +The random number generator "Mersenne Twister", which is used by GLib, +is developed and originally coded by: +Makoto Matsumoto +Takuji Nishimura diff --git a/ChangeLog b/ChangeLog index 87e7995..90d5314 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,11 @@ +1999-04-09 Sebastian Wilhelmi + + * grand.c, tests/rand-test.c: New files to implement the Mersenne + Twister Pseudo Random Number Generator. + + * glib.h, AUTHORS, Makefile.am, tests/Makefile.am: Changed + accordingly. + Thu Apr 8 21:12:30 CDT 1999 Shawn T. Amundson * Released GLib 1.3.0 diff --git a/ChangeLog.pre-2-0 b/ChangeLog.pre-2-0 index 87e7995..90d5314 100644 --- a/ChangeLog.pre-2-0 +++ b/ChangeLog.pre-2-0 @@ -1,3 +1,11 @@ +1999-04-09 Sebastian Wilhelmi + + * grand.c, tests/rand-test.c: New files to implement the Mersenne + Twister Pseudo Random Number Generator. + + * glib.h, AUTHORS, Makefile.am, tests/Makefile.am: Changed + accordingly. + Thu Apr 8 21:12:30 CDT 1999 Shawn T. Amundson * Released GLib 1.3.0 diff --git a/ChangeLog.pre-2-10 b/ChangeLog.pre-2-10 index 87e7995..90d5314 100644 --- a/ChangeLog.pre-2-10 +++ b/ChangeLog.pre-2-10 @@ -1,3 +1,11 @@ +1999-04-09 Sebastian Wilhelmi + + * grand.c, tests/rand-test.c: New files to implement the Mersenne + Twister Pseudo Random Number Generator. + + * glib.h, AUTHORS, Makefile.am, tests/Makefile.am: Changed + accordingly. + Thu Apr 8 21:12:30 CDT 1999 Shawn T. Amundson * Released GLib 1.3.0 diff --git a/ChangeLog.pre-2-12 b/ChangeLog.pre-2-12 index 87e7995..90d5314 100644 --- a/ChangeLog.pre-2-12 +++ b/ChangeLog.pre-2-12 @@ -1,3 +1,11 @@ +1999-04-09 Sebastian Wilhelmi + + * grand.c, tests/rand-test.c: New files to implement the Mersenne + Twister Pseudo Random Number Generator. + + * glib.h, AUTHORS, Makefile.am, tests/Makefile.am: Changed + accordingly. + Thu Apr 8 21:12:30 CDT 1999 Shawn T. Amundson * Released GLib 1.3.0 diff --git a/ChangeLog.pre-2-2 b/ChangeLog.pre-2-2 index 87e7995..90d5314 100644 --- a/ChangeLog.pre-2-2 +++ b/ChangeLog.pre-2-2 @@ -1,3 +1,11 @@ +1999-04-09 Sebastian Wilhelmi + + * grand.c, tests/rand-test.c: New files to implement the Mersenne + Twister Pseudo Random Number Generator. + + * glib.h, AUTHORS, Makefile.am, tests/Makefile.am: Changed + accordingly. + Thu Apr 8 21:12:30 CDT 1999 Shawn T. Amundson * Released GLib 1.3.0 diff --git a/ChangeLog.pre-2-4 b/ChangeLog.pre-2-4 index 87e7995..90d5314 100644 --- a/ChangeLog.pre-2-4 +++ b/ChangeLog.pre-2-4 @@ -1,3 +1,11 @@ +1999-04-09 Sebastian Wilhelmi + + * grand.c, tests/rand-test.c: New files to implement the Mersenne + Twister Pseudo Random Number Generator. + + * glib.h, AUTHORS, Makefile.am, tests/Makefile.am: Changed + accordingly. + Thu Apr 8 21:12:30 CDT 1999 Shawn T. Amundson * Released GLib 1.3.0 diff --git a/ChangeLog.pre-2-6 b/ChangeLog.pre-2-6 index 87e7995..90d5314 100644 --- a/ChangeLog.pre-2-6 +++ b/ChangeLog.pre-2-6 @@ -1,3 +1,11 @@ +1999-04-09 Sebastian Wilhelmi + + * grand.c, tests/rand-test.c: New files to implement the Mersenne + Twister Pseudo Random Number Generator. + + * glib.h, AUTHORS, Makefile.am, tests/Makefile.am: Changed + accordingly. + Thu Apr 8 21:12:30 CDT 1999 Shawn T. Amundson * Released GLib 1.3.0 diff --git a/ChangeLog.pre-2-8 b/ChangeLog.pre-2-8 index 87e7995..90d5314 100644 --- a/ChangeLog.pre-2-8 +++ b/ChangeLog.pre-2-8 @@ -1,3 +1,11 @@ +1999-04-09 Sebastian Wilhelmi + + * grand.c, tests/rand-test.c: New files to implement the Mersenne + Twister Pseudo Random Number Generator. + + * glib.h, AUTHORS, Makefile.am, tests/Makefile.am: Changed + accordingly. + Thu Apr 8 21:12:30 CDT 1999 Shawn T. Amundson * Released GLib 1.3.0 diff --git a/Makefile.am b/Makefile.am index bbcb360..32014e4 100644 --- a/Makefile.am +++ b/Makefile.am @@ -44,6 +44,7 @@ libglib_la_SOURCES = \ gprimes.c \ gqueue.c \ grel.c \ + grand.c \ gscanner.c \ gslist.c \ gstack.c \ diff --git a/glib.h b/glib.h index 01736c2..a772fb5 100644 --- a/glib.h +++ b/glib.h @@ -2369,6 +2369,49 @@ gpointer g_tuples_index (GTuples *tuples, gint field); +/* GRand - a good and fast random number generator: Mersenne Twister + * see http://www.math.keio.ac.jp/~matumoto/emt.html for more info. + * The range functions return a value in the intervall [min,max). + * int -> [0..2^32-1] + * int_range -> [min..max-1] + * double -> [0..1) + * double_range -> [min..max) + */ + +typedef struct _GRand GRand; +GRand* g_rand_new_with_seed (guint32 seed); +GRand* g_rand_new (); +void g_rand_free (GRand *rand); + +void g_rand_set_seed (GRand *rand, + guint32 seed); +guint32 g_rand_int (GRand *rand); +gint32 g_rand_int_range (GRand *rand, + gint32 min, + gint32 max); +gdouble g_rand_double (GRand *rand); +gdouble g_rand_double_range (GRand *rand, + gdouble min, + gdouble max); +/* This might go in, if -lm is no problem for you guys +gdouble g_rand_normal (GRand *rand, + gdouble mean, + gdouble standard_deviation); +*/ + +void g_random_set_seed (guint32 seed); +guint32 g_random_int (); +gint32 g_random_int_range (gint32 min, + gint32 max); +gdouble g_random_double (); +gdouble g_random_double_range (gdouble min, + gdouble max); +/* dito +gdouble g_random_normal (gdouble mean, + gdouble standard_deviation); +*/ + + /* Prime numbers. */ diff --git a/glib/Makefile.am b/glib/Makefile.am index bbcb360..32014e4 100644 --- a/glib/Makefile.am +++ b/glib/Makefile.am @@ -44,6 +44,7 @@ libglib_la_SOURCES = \ gprimes.c \ gqueue.c \ grel.c \ + grand.c \ gscanner.c \ gslist.c \ gstack.c \ diff --git a/glib/glib.h b/glib/glib.h index 01736c2..a772fb5 100644 --- a/glib/glib.h +++ b/glib/glib.h @@ -2369,6 +2369,49 @@ gpointer g_tuples_index (GTuples *tuples, gint field); +/* GRand - a good and fast random number generator: Mersenne Twister + * see http://www.math.keio.ac.jp/~matumoto/emt.html for more info. + * The range functions return a value in the intervall [min,max). + * int -> [0..2^32-1] + * int_range -> [min..max-1] + * double -> [0..1) + * double_range -> [min..max) + */ + +typedef struct _GRand GRand; +GRand* g_rand_new_with_seed (guint32 seed); +GRand* g_rand_new (); +void g_rand_free (GRand *rand); + +void g_rand_set_seed (GRand *rand, + guint32 seed); +guint32 g_rand_int (GRand *rand); +gint32 g_rand_int_range (GRand *rand, + gint32 min, + gint32 max); +gdouble g_rand_double (GRand *rand); +gdouble g_rand_double_range (GRand *rand, + gdouble min, + gdouble max); +/* This might go in, if -lm is no problem for you guys +gdouble g_rand_normal (GRand *rand, + gdouble mean, + gdouble standard_deviation); +*/ + +void g_random_set_seed (guint32 seed); +guint32 g_random_int (); +gint32 g_random_int_range (gint32 min, + gint32 max); +gdouble g_random_double (); +gdouble g_random_double_range (gdouble min, + gdouble max); +/* dito +gdouble g_random_normal (gdouble mean, + gdouble standard_deviation); +*/ + + /* Prime numbers. */ diff --git a/glib/grand.c b/glib/grand.c new file mode 100644 index 0000000..2725d9f --- /dev/null +++ b/glib/grand.c @@ -0,0 +1,334 @@ +/* GLIB - Library of useful routines for C programming + * Copyright (C) 1995-1997 Peter Mattis, Spencer Kimball and Josh MacDonald + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Library General Public + * License as published by the Free Software Foundation; either + * version 2 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with this library; if not, write to the + * Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ + +/* Originally developed and coded by Makoto Matsumoto and Takuji + * Nishimura. Please mail , if you're using + * code from this file in your own programs or libraries. + * Further information on the Mersenne Twister can be found at + * http://www.math.keio.ac.jp/~matumoto/emt.html + */ + +/* + * Modified by the GLib Team and others 1997-1999. See the AUTHORS + * file for a list of people on the GLib Team. See the ChangeLog + * files for a list of changes. These files are distributed with + * GLib at ftp://ftp.gtk.org/pub/gtk/. + */ + +#include +#include +#include + +G_LOCK_DEFINE_STATIC (global_random); +static GRand* global_random = NULL; + +/* Period parameters */ +#define N 624 +#define M 397 +#define MATRIX_A 0x9908b0df /* constant vector a */ +#define UPPER_MASK 0x80000000 /* most significant w-r bits */ +#define LOWER_MASK 0x7fffffff /* least significant r bits */ + +/* Tempering parameters */ +#define TEMPERING_MASK_B 0x9d2c5680 +#define TEMPERING_MASK_C 0xefc60000 +#define TEMPERING_SHIFT_U(y) (y >> 11) +#define TEMPERING_SHIFT_S(y) (y << 7) +#define TEMPERING_SHIFT_T(y) (y << 15) +#define TEMPERING_SHIFT_L(y) (y >> 18) + +struct _GRand +{ + guint32 mt[N]; /* the array for the state vector */ + guint mti; + gboolean have_next_normal; + gdouble next_normal; +}; + +GRand* +g_rand_new_with_seed (guint32 seed) +{ + GRand *rand = g_new0 (GRand, 1); + g_rand_set_seed (rand, seed); + return rand; +} + +GRand* +g_rand_new () +{ + guint32 seed = 0; + GTimeVal now; + FILE* dev_random = fopen("/dev/random", "rb"); + + if (dev_random) + { + if (fread (&seed, sizeof (seed), 1, dev_random) != 1) + seed = 0; + fclose (dev_random); + } + + /* Using /dev/random alone makes the seed computable for the + outside. This might pose security problems somewhere. This should + yield better values */ + + g_get_current_time (&now); + seed ^= now.tv_sec ^ now.tv_usec; + + return g_rand_new_with_seed (seed); +} + +void +g_rand_free (GRand* rand) +{ + g_return_if_fail (rand); + + g_free (rand); +} + +void +g_rand_set_seed (GRand* rand, guint32 seed) +{ + g_return_if_fail (rand); + + /* setting initial seeds to mt[N] using */ + /* the generator Line 25 of Table 1 in */ + /* [KNUTH 1981, The Art of Computer Programming */ + /* Vol. 2 (2nd Ed.), pp102] */ + rand->mt[0]= seed & 0xffffffff; + for (rand->mti=1; rand->mtimti++) + rand->mt[rand->mti] = (69069 * rand->mt[rand->mti-1]) & 0xffffffff; + + rand->have_next_normal = FALSE; +} + +guint32 +g_rand_int (GRand* rand) +{ + guint32 y; + static const guint32 mag01[2]={0x0, MATRIX_A}; + /* mag01[x] = x * MATRIX_A for x=0,1 */ + + g_return_val_if_fail (rand, 0); + + if (rand->mti >= N) { /* generate N words at one time */ + int kk; + + for (kk=0;kkmt[kk]&UPPER_MASK)|(rand->mt[kk+1]&LOWER_MASK); + rand->mt[kk] = rand->mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1]; + } + for (;kkmt[kk]&UPPER_MASK)|(rand->mt[kk+1]&LOWER_MASK); + rand->mt[kk] = rand->mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1]; + } + y = (rand->mt[N-1]&UPPER_MASK)|(rand->mt[0]&LOWER_MASK); + rand->mt[N-1] = rand->mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1]; + + rand->mti = 0; + } + + y = rand->mt[rand->mti++]; + y ^= TEMPERING_SHIFT_U(y); + y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B; + y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C; + y ^= TEMPERING_SHIFT_L(y); + + return y; +} + +gint32 +g_rand_int_range (GRand* rand, gint32 min, gint32 max) +{ + guint32 dist = max - min; + guint32 random; + + g_return_val_if_fail (rand, min); + g_return_val_if_fail (max > min, min); + + if (dist <= 0x10000L) /* 2^16 */ + { + /* All tricks doing modulo calculations do not have a good + distribution -> We must use this slower method for maximal + quality, but this method is only good for (max - min) <= 2^16 */ + + random = (gint32) g_rand_double_range (rand, 0, dist); + /* we'd rather use the following, if -lm is allowed later on: + random = (gint32) floor (g_rand_double_range (rand, 0, dist)); */ + } + else + { + /* Now it's harder to make it right. We calculate the smallest m, + such that dist < 2 ^ m, then we calculate a random number in + [1..2^32-1] and rightshift it by 32 - m. Then we test, if it + is smaller than dist and if not, get a new number and so + forth until we get a number smaller than dist. We just return + this. */ + guint32 border = 0x20000L; /* 2^17 */ + guint right_shift = 15; /* 32 - 17 */ + + if (dist >= 0x80000000) /* in the case of dist > 2^31 our loop + below will be infinite */ + { + right_shift = 0; + } + else + { + while (dist >= border) + { + border <<= 1; + right_shift--; + } + } + do + { + random = g_rand_int (rand) >> right_shift; + } while (random >= dist); + } + return min + random; +} + +/* transform [0..2^32-1] -> [0..1) */ +#define G_RAND_DOUBLE_TRANSFORM 2.3283064365386963e-10 + +gdouble +g_rand_double (GRand* rand) +{ + return g_rand_int (rand) * G_RAND_DOUBLE_TRANSFORM; +} + +gdouble +g_rand_double_range (GRand* rand, gdouble min, gdouble max) +{ + return g_rand_int (rand) * ((max - min) * G_RAND_DOUBLE_TRANSFORM) + min; +} + + +#if WE_REALLY_WANT_HAVE_MATH_LIB_LINKED +gdouble +g_rand_normal (GRand* rand, gdouble mean, gdouble standard_deviation) +{ + /* For a description of the used algorithm see Knuth: "The Art of + Computer Programming", Vol.2, Second Edition, Page 117: Polar + method for normal deviates due to Box, Muller, Marsaglia */ + gdouble normal; + g_return_val_if_fail (rand, 0); + + if (rand->have_next_normal) + { + rand->have_next_normal = FALSE; + normal = rand->next_normal; + } + else + { + gdouble u1; + gdouble u2 = g_rand_double_range (rand, -1, 1); + gdouble s, f; + do + { + u1 = u2; + u2 = g_rand_double_range (rand, -1, 1); + s = u1 * u1 + u2 * u2; + } while (s >= 1.0); + f = sqrt (-2 * log (s) / s); + normal = u1 * f; + rand->next_normal = u2 * f; + rand->have_next_normal = TRUE; + } + return mean + normal * standard_deviation; +} +#endif + +guint32 +g_random_int (void) +{ + guint32 result; + G_LOCK (global_random); + if (!global_random) + global_random = g_rand_new (); + + result = g_rand_int (global_random); + G_UNLOCK (global_random); + return result; +} + +gint32 +g_random_int_range (gint32 min, gint32 max) +{ + gint32 result; + G_LOCK (global_random); + if (!global_random) + global_random = g_rand_new (); + + result = g_rand_int_range (global_random, min, max); + G_UNLOCK (global_random); + return result; +} + +gdouble +g_random_double (void) +{ + double result; + G_LOCK (global_random); + if (!global_random) + global_random = g_rand_new (); + + result = g_rand_double (global_random); + G_UNLOCK (global_random); + return result; +} + +gdouble +g_random_double_range (gdouble min, gdouble max) +{ + double result; + G_LOCK (global_random); + if (!global_random) + global_random = g_rand_new (); + + result = g_rand_double_range (global_random, min, max); + G_UNLOCK (global_random); + return result; +} + +#if WE_REALLY_WANT_HAVE_MATH_LIB_LINKED +gdouble +g_random_normal (gdouble mean, gdouble standard_deviation) +{ + double result; + G_LOCK (global_random); + if (!global_random) + global_random = g_rand_new (); + + result = g_rand_normal (global_random, mean, standard_deviation); + G_UNLOCK (global_random); + return result; +} +#endif + +void +g_random_set_seed (guint32 seed) +{ + G_LOCK (global_random); + if (!global_random) + global_random = g_rand_new_with_seed (seed); + else + g_rand_set_seed (global_random, seed); + G_UNLOCK (global_random); +} + diff --git a/grand.c b/grand.c new file mode 100644 index 0000000..2725d9f --- /dev/null +++ b/grand.c @@ -0,0 +1,334 @@ +/* GLIB - Library of useful routines for C programming + * Copyright (C) 1995-1997 Peter Mattis, Spencer Kimball and Josh MacDonald + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Library General Public + * License as published by the Free Software Foundation; either + * version 2 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with this library; if not, write to the + * Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ + +/* Originally developed and coded by Makoto Matsumoto and Takuji + * Nishimura. Please mail , if you're using + * code from this file in your own programs or libraries. + * Further information on the Mersenne Twister can be found at + * http://www.math.keio.ac.jp/~matumoto/emt.html + */ + +/* + * Modified by the GLib Team and others 1997-1999. See the AUTHORS + * file for a list of people on the GLib Team. See the ChangeLog + * files for a list of changes. These files are distributed with + * GLib at ftp://ftp.gtk.org/pub/gtk/. + */ + +#include +#include +#include + +G_LOCK_DEFINE_STATIC (global_random); +static GRand* global_random = NULL; + +/* Period parameters */ +#define N 624 +#define M 397 +#define MATRIX_A 0x9908b0df /* constant vector a */ +#define UPPER_MASK 0x80000000 /* most significant w-r bits */ +#define LOWER_MASK 0x7fffffff /* least significant r bits */ + +/* Tempering parameters */ +#define TEMPERING_MASK_B 0x9d2c5680 +#define TEMPERING_MASK_C 0xefc60000 +#define TEMPERING_SHIFT_U(y) (y >> 11) +#define TEMPERING_SHIFT_S(y) (y << 7) +#define TEMPERING_SHIFT_T(y) (y << 15) +#define TEMPERING_SHIFT_L(y) (y >> 18) + +struct _GRand +{ + guint32 mt[N]; /* the array for the state vector */ + guint mti; + gboolean have_next_normal; + gdouble next_normal; +}; + +GRand* +g_rand_new_with_seed (guint32 seed) +{ + GRand *rand = g_new0 (GRand, 1); + g_rand_set_seed (rand, seed); + return rand; +} + +GRand* +g_rand_new () +{ + guint32 seed = 0; + GTimeVal now; + FILE* dev_random = fopen("/dev/random", "rb"); + + if (dev_random) + { + if (fread (&seed, sizeof (seed), 1, dev_random) != 1) + seed = 0; + fclose (dev_random); + } + + /* Using /dev/random alone makes the seed computable for the + outside. This might pose security problems somewhere. This should + yield better values */ + + g_get_current_time (&now); + seed ^= now.tv_sec ^ now.tv_usec; + + return g_rand_new_with_seed (seed); +} + +void +g_rand_free (GRand* rand) +{ + g_return_if_fail (rand); + + g_free (rand); +} + +void +g_rand_set_seed (GRand* rand, guint32 seed) +{ + g_return_if_fail (rand); + + /* setting initial seeds to mt[N] using */ + /* the generator Line 25 of Table 1 in */ + /* [KNUTH 1981, The Art of Computer Programming */ + /* Vol. 2 (2nd Ed.), pp102] */ + rand->mt[0]= seed & 0xffffffff; + for (rand->mti=1; rand->mtimti++) + rand->mt[rand->mti] = (69069 * rand->mt[rand->mti-1]) & 0xffffffff; + + rand->have_next_normal = FALSE; +} + +guint32 +g_rand_int (GRand* rand) +{ + guint32 y; + static const guint32 mag01[2]={0x0, MATRIX_A}; + /* mag01[x] = x * MATRIX_A for x=0,1 */ + + g_return_val_if_fail (rand, 0); + + if (rand->mti >= N) { /* generate N words at one time */ + int kk; + + for (kk=0;kkmt[kk]&UPPER_MASK)|(rand->mt[kk+1]&LOWER_MASK); + rand->mt[kk] = rand->mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1]; + } + for (;kkmt[kk]&UPPER_MASK)|(rand->mt[kk+1]&LOWER_MASK); + rand->mt[kk] = rand->mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1]; + } + y = (rand->mt[N-1]&UPPER_MASK)|(rand->mt[0]&LOWER_MASK); + rand->mt[N-1] = rand->mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1]; + + rand->mti = 0; + } + + y = rand->mt[rand->mti++]; + y ^= TEMPERING_SHIFT_U(y); + y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B; + y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C; + y ^= TEMPERING_SHIFT_L(y); + + return y; +} + +gint32 +g_rand_int_range (GRand* rand, gint32 min, gint32 max) +{ + guint32 dist = max - min; + guint32 random; + + g_return_val_if_fail (rand, min); + g_return_val_if_fail (max > min, min); + + if (dist <= 0x10000L) /* 2^16 */ + { + /* All tricks doing modulo calculations do not have a good + distribution -> We must use this slower method for maximal + quality, but this method is only good for (max - min) <= 2^16 */ + + random = (gint32) g_rand_double_range (rand, 0, dist); + /* we'd rather use the following, if -lm is allowed later on: + random = (gint32) floor (g_rand_double_range (rand, 0, dist)); */ + } + else + { + /* Now it's harder to make it right. We calculate the smallest m, + such that dist < 2 ^ m, then we calculate a random number in + [1..2^32-1] and rightshift it by 32 - m. Then we test, if it + is smaller than dist and if not, get a new number and so + forth until we get a number smaller than dist. We just return + this. */ + guint32 border = 0x20000L; /* 2^17 */ + guint right_shift = 15; /* 32 - 17 */ + + if (dist >= 0x80000000) /* in the case of dist > 2^31 our loop + below will be infinite */ + { + right_shift = 0; + } + else + { + while (dist >= border) + { + border <<= 1; + right_shift--; + } + } + do + { + random = g_rand_int (rand) >> right_shift; + } while (random >= dist); + } + return min + random; +} + +/* transform [0..2^32-1] -> [0..1) */ +#define G_RAND_DOUBLE_TRANSFORM 2.3283064365386963e-10 + +gdouble +g_rand_double (GRand* rand) +{ + return g_rand_int (rand) * G_RAND_DOUBLE_TRANSFORM; +} + +gdouble +g_rand_double_range (GRand* rand, gdouble min, gdouble max) +{ + return g_rand_int (rand) * ((max - min) * G_RAND_DOUBLE_TRANSFORM) + min; +} + + +#if WE_REALLY_WANT_HAVE_MATH_LIB_LINKED +gdouble +g_rand_normal (GRand* rand, gdouble mean, gdouble standard_deviation) +{ + /* For a description of the used algorithm see Knuth: "The Art of + Computer Programming", Vol.2, Second Edition, Page 117: Polar + method for normal deviates due to Box, Muller, Marsaglia */ + gdouble normal; + g_return_val_if_fail (rand, 0); + + if (rand->have_next_normal) + { + rand->have_next_normal = FALSE; + normal = rand->next_normal; + } + else + { + gdouble u1; + gdouble u2 = g_rand_double_range (rand, -1, 1); + gdouble s, f; + do + { + u1 = u2; + u2 = g_rand_double_range (rand, -1, 1); + s = u1 * u1 + u2 * u2; + } while (s >= 1.0); + f = sqrt (-2 * log (s) / s); + normal = u1 * f; + rand->next_normal = u2 * f; + rand->have_next_normal = TRUE; + } + return mean + normal * standard_deviation; +} +#endif + +guint32 +g_random_int (void) +{ + guint32 result; + G_LOCK (global_random); + if (!global_random) + global_random = g_rand_new (); + + result = g_rand_int (global_random); + G_UNLOCK (global_random); + return result; +} + +gint32 +g_random_int_range (gint32 min, gint32 max) +{ + gint32 result; + G_LOCK (global_random); + if (!global_random) + global_random = g_rand_new (); + + result = g_rand_int_range (global_random, min, max); + G_UNLOCK (global_random); + return result; +} + +gdouble +g_random_double (void) +{ + double result; + G_LOCK (global_random); + if (!global_random) + global_random = g_rand_new (); + + result = g_rand_double (global_random); + G_UNLOCK (global_random); + return result; +} + +gdouble +g_random_double_range (gdouble min, gdouble max) +{ + double result; + G_LOCK (global_random); + if (!global_random) + global_random = g_rand_new (); + + result = g_rand_double_range (global_random, min, max); + G_UNLOCK (global_random); + return result; +} + +#if WE_REALLY_WANT_HAVE_MATH_LIB_LINKED +gdouble +g_random_normal (gdouble mean, gdouble standard_deviation) +{ + double result; + G_LOCK (global_random); + if (!global_random) + global_random = g_rand_new (); + + result = g_rand_normal (global_random, mean, standard_deviation); + G_UNLOCK (global_random); + return result; +} +#endif + +void +g_random_set_seed (guint32 seed) +{ + G_LOCK (global_random); + if (!global_random) + global_random = g_rand_new_with_seed (seed); + else + g_rand_set_seed (global_random, seed); + G_UNLOCK (global_random); +} + diff --git a/gthread/Makefile.am b/gthread/Makefile.am index 2752455..f24dffb 100644 --- a/gthread/Makefile.am +++ b/gthread/Makefile.am @@ -22,4 +22,4 @@ libgthread_la_LDFLAGS = \ libgthread_la_LIBADD = @G_THREAD_LIBS@ noinst_PROGRAMS = testgthread -testgthread_LDADD = ../libglib.la libgthread.la +testgthread_LDADD = ../libglib.la libgthread.la diff --git a/tests/Makefile.am b/tests/Makefile.am index 62fcec0..539d2d1 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -9,6 +9,7 @@ TESTS = \ list-test \ node-test \ queue-test \ + rand-test \ relation-test \ slist-test \ stack-test \ @@ -26,6 +27,7 @@ hash_test_LDADD = $(top_builddir)/libglib.la list_test_LDADD = $(top_builddir)/libglib.la node_test_LDADD = $(top_builddir)/libglib.la queue_test_LDADD = $(top_builddir)/libglib.la +rand_test_LDADD = $(top_builddir)/libglib.la relation_test_LDADD = $(top_builddir)/libglib.la slist_test_LDADD = $(top_builddir)/libglib.la stack_test_LDADD = $(top_builddir)/libglib.la diff --git a/tests/rand-test.c b/tests/rand-test.c new file mode 100644 index 0000000..d3af9eb --- /dev/null +++ b/tests/rand-test.c @@ -0,0 +1,43 @@ +#include + +const gint32 first_numbers[] = +{ + 0x7a7a7a7a, + 0x20aea82a, + 0xcab337ab, + 0xdcf770ea, + 0xdf552b2f, + 0x32d1ef7f, + 0x6bed6dd9, + 0x7222df44, + 0x6b842128, + 0x07f8579a, + 0x9dad1004, + 0x2df264f2, + 0x13b48989, + 0xf2929475, + 0x30f30c97, + 0x3f9a1ea7, + 0x3bf04710, + 0xb85bd69e, + 0x790a48b0, + 0xfa06b85f, + 0xa64cc9e3 +}; + +const gint length = sizeof (first_numbers) / sizeof (first_numbers[0]); + +int main() +{ + guint i; + + GRand* rand = g_rand_new_with_seed (first_numbers[0]); + + for (i = 1; i < length; i++) + g_assert (first_numbers[i]); + + g_rand_free (rand); + + return 0; +} + -- 2.7.4