diff options
author | Søren Sandmann Pedersen <sandmann@daimi.au.dk> | 2009-08-30 09:03:26 -0400 |
---|---|---|
committer | Søren Sandmann Pedersen <sandmann@daimi.au.dk> | 2009-08-30 09:03:26 -0400 |
commit | d9b2da69993c66b8d7039c15cbe87539e955a6d8 (patch) | |
tree | 7d6907c978094729d4cfb7eff37dabb45ecef7db | |
parent | 8e470dd69047ed08963b8557fb6251b0cbd47b52 (diff) |
Add a random number generator from glib
-rw-r--r-- | Makefile.am | 1 | ||||
-rw-r--r-- | libnul.h | 14 | ||||
-rw-r--r-- | random.c | 502 |
3 files changed, 517 insertions, 0 deletions
diff --git a/Makefile.am b/Makefile.am index abc5e53..e999729 100644 --- a/Makefile.am +++ b/Makefile.am @@ -10,6 +10,7 @@ libnul_la_SOURCES = \ array.c \ epoll.c \ dbus-int.c \ + random.c \ signal-handler.c \ file-utils.c \ libnul.h @@ -307,6 +307,20 @@ nul_bool_t nul_fd_is_watched (int fd); */ char *nul_canonicalize_filename (const char *filename); + +/* + * Random numbers + */ +void nul_random_set_seed (uint32_t seed); +uint32_t nul_random_int (void); +int32_t nul_random_int_range (int32_t begin, + int32_t end); +double nul_random_double (void); +double nul_random_double_range (double begin, + double end); + +#define nul_random_boolean() ((nul_random_int () & (1 << 15)) != 0) + /* * Unix signal handlers */ diff --git a/random.c b/random.c new file mode 100644 index 0000000..cc84023 --- /dev/null +++ b/random.c @@ -0,0 +1,502 @@ +/* libnul - another utility library */ + +/* Based on grand.c from GLib. GLib copyright: + * + * 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 Lesser 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 + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser 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 <matumoto@math.keio.ac.jp>, 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.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html + * This code was adapted to glib by Sebastian Wilhelmi. + */ + +/* + * Modified by the GLib Team and others 1997-2000. 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/. + */ + +/* + * MT safe + */ + +#include "config.h" +#include "libnul.h" + +#include <math.h> +#include <errno.h> +#include <stdio.h> +#include <string.h> +#include <sys/types.h> +#ifdef HAVE_UNISTD_H +#include <unistd.h> +#endif + +typedef struct _rand_t rand_t; + +G_LOCK_DEFINE_STATIC (global_random); +static rand_t* 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 _rand_t +{ + uint32_t mt[N]; /* the array for the state vector */ + guint mti; +}; + +/** + * nul_rand_set_seed: + * @rand_: a #rand_t. + * @seed: a value to reinitialize the random number generator. + * + * Sets the seed for the random number generator #rand_t to @seed. + **/ +static void +nul_rand_set_seed (rand_t* rand, uint32_t seed) +{ + g_return_if_fail (rand != NULL); + + /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ + /* In the previous version (see above), MSBs of the */ + /* seed affect only MSBs of the array mt[]. */ + + rand->mt[0]= seed; + for (rand->mti=1; rand->mti<N; rand->mti++) + rand->mt[rand->mti] = 1812433253UL * + (rand->mt[rand->mti-1] ^ (rand->mt[rand->mti-1] >> 30)) + rand->mti; +} + +/** + * nul_rand_set_seed_array: + * @rand_: a #rand_t. + * @seed: array to initialize with + * @seed_length: length of array + * + * Initializes the random number generator by an array of + * longs. Array can be of arbitrary size, though only the + * first 624 values are taken. This function is useful + * if you have many low entropy seeds, or if you require more then + * 32bits of actual entropy for your application. + * + * Since: 2.4 + **/ +static void +nul_rand_set_seed_array (rand_t* rand, const uint32_t *seed, guint seed_length) +{ + int i, j, k; + + g_return_if_fail (rand != NULL); + g_return_if_fail (seed_length >= 1); + + nul_rand_set_seed (rand, 19650218UL); + + i=1; j=0; + k = (N>seed_length ? N : seed_length); + for (; k; k--) + { + rand->mt[i] = (rand->mt[i] ^ + ((rand->mt[i-1] ^ (rand->mt[i-1] >> 30)) * 1664525UL)) + + seed[j] + j; /* non linear */ + rand->mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ + i++; j++; + if (i>=N) + { + rand->mt[0] = rand->mt[N-1]; + i=1; + } + if (j>=seed_length) + j=0; + } + for (k=N-1; k; k--) + { + rand->mt[i] = (rand->mt[i] ^ + ((rand->mt[i-1] ^ (rand->mt[i-1] >> 30)) * 1566083941UL)) + - i; /* non linear */ + rand->mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ + i++; + if (i>=N) + { + rand->mt[0] = rand->mt[N-1]; + i=1; + } + } + + rand->mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ +} + +/** + * nul_rand_new_with_seed: + * @seed: a value to initialize the random number generator. + * + * Creates a new random number generator initialized with @seed. + * + * Return value: the new #rand_t. + **/ +static rand_t* +nul_rand_new_with_seed (uint32_t seed) +{ + rand_t *rand = g_new0 (rand_t, 1); + nul_rand_set_seed (rand, seed); + return rand; +} + +/** + * nul_rand_new_with_seed_array: + * @seed: an array of seeds to initialize the random number generator. + * @seed_length: an array of seeds to initialize the random number generator. + * + * Creates a new random number generator initialized with @seed. + * + * Return value: the new #rand_t. + * + * Since: 2.4 + **/ +static rand_t* +nul_rand_new_with_seed_array (const uint32_t *seed, guint seed_length) +{ + rand_t *rand = g_new0 (rand_t, 1); + nul_rand_set_seed_array (rand, seed, seed_length); + return rand; +} + +/** + * nul_rand_new: + * + * Creates a new random number generator initialized with a seed taken + * either from <filename>/dev/urandom</filename> (if existing) or from + * the current time (as a fallback). + * + * Return value: the new #rand_t. + **/ +static rand_t* +nul_rand_new (void) +{ + uint32_t seed[4]; + GTimeVal now; +#ifdef G_OS_UNIX + static gboolean dev_urandom_exists = TRUE; + + if (dev_urandom_exists) + { + FILE* dev_urandom; + + do + { + errno = 0; + dev_urandom = fopen("/dev/urandom", "rb"); + } + while G_UNLIKELY (errno == EINTR); + + if (dev_urandom) + { + int r; + + do + { + errno = 0; + r = fread (seed, sizeof (seed), 1, dev_urandom); + } + while G_UNLIKELY (errno == EINTR); + + if (r != 1) + dev_urandom_exists = FALSE; + + fclose (dev_urandom); + } + else + dev_urandom_exists = FALSE; + } +#else + static gboolean dev_urandom_exists = FALSE; +#endif + + if (!dev_urandom_exists) + { + g_get_current_time (&now); + seed[0] = now.tv_sec; + seed[1] = now.tv_usec; + seed[2] = getpid (); +#ifdef G_OS_UNIX + seed[3] = getppid (); +#else + seed[3] = 0; +#endif + } + + return nul_rand_new_with_seed_array (seed, 4); +} + +/** + * nul_rand_int: + * @rand_: a #rand_t. + * + * Returns the next random #uint32_t from @rand_ equally distributed over + * the range [0..2^32-1]. + * + * Return value: A random number. + **/ +static uint32_t +nul_rand_int (rand_t* rand) +{ + uint32_t y; + static const uint32_t mag01[2]={0x0, MATRIX_A}; + /* mag01[x] = x * MATRIX_A for x=0,1 */ + + g_return_val_if_fail (rand != NULL, 0); + + if (rand->mti >= N) { /* generate N words at one time */ + int kk; + + for (kk=0;kk<N-M;kk++) { + y = (rand->mt[kk]&UPPER_MASK)|(rand->mt[kk+1]&LOWER_MASK); + rand->mt[kk] = rand->mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1]; + } + for (;kk<N-1;kk++) { + y = (rand->mt[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; +} + +/* transform [0..2^32] -> [0..1] */ +#define NUL_RAND_DOUBLE_TRANSFORM 2.3283064365386962890625e-10 + +/** + * nul_rand_int_range: + * @rand_: a #rand_t. + * @begin: lower closed bound of the interval. + * @end: upper open bound of the interval. + * + * Returns the next random #int32_t from @rand_ equally distributed over + * the range [@begin..@end-1]. + * + * Return value: A random number. + **/ +static int32_t +nul_rand_int_range (rand_t* rand, int32_t begin, int32_t end) +{ + uint32_t dist = end - begin; + uint32_t random; + + g_return_val_if_fail (rand != NULL, begin); + g_return_val_if_fail (end > begin, begin); + + if (dist == 0) + random = 0; + else + { + /* maxvalue is set to the predecessor of the greatest + * multiple of dist less or equal 2^32. */ + uint32_t maxvalue; + if (dist <= 0x80000000u) /* 2^31 */ + { + /* maxvalue = 2^32 - 1 - (2^32 % dist) */ + uint32_t leftover = (0x80000000u % dist) * 2; + if (leftover >= dist) leftover -= dist; + maxvalue = 0xffffffffu - leftover; + } + else + maxvalue = dist - 1; + + do + random = nul_rand_int (rand); + while (random > maxvalue); + + random %= dist; + } + return begin + random; +} + +/** + * nul_rand_double: + * @rand_: a #rand_t. + * + * Returns the next random #double from @rand_ equally distributed over + * the range [0..1). + * + * Return value: A random number. + **/ +static double +nul_rand_double (rand_t* rand) +{ + /* We set all 52 bits after the point for this, not only the first + 32. Thats why we need two calls to nul_rand_int */ + double retval = nul_rand_int (rand) * NUL_RAND_DOUBLE_TRANSFORM; + retval = (retval + nul_rand_int (rand)) * NUL_RAND_DOUBLE_TRANSFORM; + + /* The following might happen due to very bad rounding luck, but + * actually this should be more than rare, we just try again then */ + if (retval >= 1.0) + return nul_rand_double (rand); + + return retval; +} + +/** + * nul_rand_double_range: + * @rand_: a #rand_t. + * @begin: lower closed bound of the interval. + * @end: upper open bound of the interval. + * + * Returns the next random #double from @rand_ equally distributed over + * the range [@begin..@end). + * + * Return value: A random number. + **/ +static double +nul_rand_double_range (rand_t* rand, double begin, double end) +{ + return nul_rand_double (rand) * (end - begin) + begin; +} + +/** + * nul_random_int: + * + * Return a random #uint32_t equally distributed over the range + * [0..2^32-1]. + * + * Return value: A random number. + **/ +uint32_t +nul_random_int (void) +{ + uint32_t result; + G_LOCK (global_random); + if (!global_random) + global_random = nul_rand_new (); + + result = nul_rand_int (global_random); + G_UNLOCK (global_random); + return result; +} + +/** + * nul_random_int_range: + * @begin: lower closed bound of the interval. + * @end: upper open bound of the interval. + * + * Returns a random #int32_t equally distributed over the range + * [@begin..@end-1]. + * + * Return value: A random number. + **/ +int32_t +nul_random_int_range (int32_t begin, int32_t end) +{ + int32_t result; + G_LOCK (global_random); + if (!global_random) + global_random = nul_rand_new (); + + result = nul_rand_int_range (global_random, begin, end); + G_UNLOCK (global_random); + return result; +} + +/** + * nul_random_double: + * + * Returns a random #double equally distributed over the range [0..1). + * + * Return value: A random number. + **/ +double +nul_random_double (void) +{ + double result; + G_LOCK (global_random); + if (!global_random) + global_random = nul_rand_new (); + + result = nul_rand_double (global_random); + G_UNLOCK (global_random); + return result; +} + +/** + * nul_random_double_range: + * @begin: lower closed bound of the interval. + * @end: upper open bound of the interval. + * + * Returns a random #double equally distributed over the range [@begin..@end). + * + * Return value: A random number. + **/ +double +nul_random_double_range (double begin, double end) +{ + double result; + G_LOCK (global_random); + if (!global_random) + global_random = nul_rand_new (); + + result = nul_rand_double_range (global_random, begin, end); + G_UNLOCK (global_random); + return result; +} + +/** + * nul_random_set_seed: + * @seed: a value to reinitialize the global random number generator. + * + * Sets the seed for the global random number generator, which is used + * by the <function>nul_random_*</function> functions, to @seed. + **/ +void +nul_random_set_seed (uint32_t seed) +{ + G_LOCK (global_random); + if (!global_random) + global_random = nul_rand_new_with_seed (seed); + else + nul_rand_set_seed (global_random, seed); + G_UNLOCK (global_random); +} + + +#define __NUL_RAND_C__ |