grand.c 9.92 KB
Newer Older
1 2 3 4
/* 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
5
 * modify it under the terms of the GNU Lesser General Public
6 7 8 9 10 11
 * 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
12
 * Lesser General Public License for more details.
13
 *
14
 * You should have received a copy of the GNU Lesser General Public
15 16 17 18 19 20 21 22 23 24
 * 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.keio.ac.jp/~matumoto/emt.html
25
 * This code was adapted to glib by Sebastian Wilhelmi <wilhelmi@ira.uka.de>.
26 27 28
 */

/*
29
 * Modified by the GLib Team and others 1997-2000.  See the AUTHORS
30 31 32 33 34
 * 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/.  
 */

35 36 37 38
/* 
 * MT safe
 */

39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
#include <glib.h>
#include <math.h>
#include <stdio.h>

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; 
};

67 68 69 70 71 72 73 74
/**
 * g_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 #GRand.
 **/
75 76 77 78 79 80 81 82
GRand*
g_rand_new_with_seed (guint32 seed)
{
  GRand *rand = g_new0 (GRand, 1);
  g_rand_set_seed (rand, seed);
  return rand;
}

83 84 85 86 87 88 89 90 91
/**
 * g_rand_new:
 * 
 * Creates a new random number generator initialized with a seed taken
 * either from /dev/urandom (if existing) or from the current time (as
 * a fallback).
 * 
 * Return value: the new #GRand.
 **/
92
GRand* 
93
g_rand_new (void)
94
{
95
  guint32 seed;
96
  GTimeVal now;
97
  static gboolean dev_urandom_exists = TRUE;
98
  
99
  if (dev_urandom_exists)
100
    {
101 102
      FILE* dev_urandom = fopen("/dev/urandom", "rb");
      if (dev_urandom)
103
	{
104 105 106
	  if (fread (&seed, sizeof (seed), 1, dev_urandom) != 1)
	    dev_urandom_exists = FALSE;
	  fclose (dev_urandom);
107 108
	}	
      else
109 110 111 112 113 114
	dev_urandom_exists = FALSE;
    }
  if (!dev_urandom_exists)
    {  
      g_get_current_time (&now);
      seed = now.tv_sec ^ now.tv_usec;
115 116 117 118 119
    }

  return g_rand_new_with_seed (seed);
}

120 121 122 123 124 125
/**
 * g_rand_free:
 * @rand: a #GRand.
 *
 * Frees the memory allocated for the #GRand.
 **/
126 127 128
void
g_rand_free (GRand* rand)
{
129
  g_return_if_fail (rand != NULL);
130 131 132 133

  g_free (rand);
}

134 135 136 137 138 139 140
/**
 * g_rand_set_seed:
 * @rand: a #GRand.
 * @seed: a value to reinitialize the random number generator.
 *
 * Sets the seed for the random number generator #GRand to @seed.
 **/
141 142 143
void
g_rand_set_seed (GRand* rand, guint32 seed)
{
144
  g_return_if_fail (rand != NULL);
145 146 147 148 149

  /* 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]                  */
150 151 152 153
  
  if (seed == 0) /* This would make the PRNG procude only zeros */
    seed = 0x6b842128; /* Just set it to another number */

154 155 156 157 158
  rand->mt[0]= seed & 0xffffffff;
  for (rand->mti=1; rand->mti<N; rand->mti++)
    rand->mt[rand->mti] = (69069 * rand->mt[rand->mti-1]) & 0xffffffff;
}

159 160 161 162 163 164 165 166 167
/**
 * g_rand_int:
 * @rand: a #GRand.
 *
 * Return the next random #guint32 from @rand equaly distributed over
 * the range [0..2^32-1].
 *
 * Return value: A random number.
 **/
168 169 170 171 172 173 174
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 */

175
  g_return_val_if_fail (rand != NULL, 0);
176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202

  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; 
}

203 204 205 206 207 208 209 210 211 212 213
/**
 * g_rand_int_range:
 * @rand: a #GRand.
 * @min: lower closed bound of the interval.
 * @max: upper open bound of the interval.
 *
 * Return the next random #gint32 from @rand equaly distributed over
 * the range [@min..@max-1].
 *
 * Return value: A random number.
 **/
214 215 216 217 218 219
gint32 
g_rand_int_range (GRand* rand, gint32 min, gint32 max)
{
  guint32 dist = max - min;
  guint32 random;

220
  g_return_val_if_fail (rand != NULL, min);
221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267
  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

268 269 270 271 272 273 274 275 276
/**
 * g_rand_double:
 * @rand: a #GRand.
 *
 * Return the next random #gdouble from @rand equaly distributed over
 * the range [0..1).
 *
 * Return value: A random number.
 **/
277 278 279 280 281 282
gdouble 
g_rand_double (GRand* rand)
{                            
  return g_rand_int (rand) * G_RAND_DOUBLE_TRANSFORM;
}

283 284 285 286 287 288 289 290 291 292 293
/**
 * g_rand_double_range:
 * @rand: a #GRand.
 * @min: lower closed bound of the interval.
 * @max: upper open bound of the interval.
 *
 * Return the next random #gdouble from @rand equaly distributed over
 * the range [@min..@max).
 *
 * Return value: A random number.
 **/
294 295 296 297 298 299
gdouble 
g_rand_double_range (GRand* rand, gdouble min, gdouble max)
{
  return g_rand_int (rand) * ((max - min) * G_RAND_DOUBLE_TRANSFORM)  + min;
}

300 301 302 303 304 305 306 307
/**
 * g_random_int:
 *
 * Return a random #guint32 equaly distributed over the range
 * [0..2^32-1].
 *
 * Return value: A random number.
 **/
308 309 310 311 312 313 314 315 316 317 318 319 320
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;
}

321 322 323 324 325 326 327 328 329 330
/**
 * g_random_int_range:
 * @min: lower closed bound of the interval.
 * @max: upper open bound of the interval.
 *
 * Return a random #gint32 equaly distributed over the range
 * [@min..@max-1].
 *
 * Return value: A random number.
 **/
331 332 333 334 335 336 337 338 339 340 341 342 343
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;
}

344 345 346 347 348 349 350
/**
 * g_random_double:
 *
 * Return a random #gdouble equaly distributed over the range [0..1).
 *
 * Return value: A random number.
 **/
351 352 353 354 355 356 357 358 359 360 361 362 363
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;
}

364 365 366 367 368 369 370 371 372
/**
 * g_random_double_range:
 * @min: lower closed bound of the interval.
 * @max: upper open bound of the interval.
 *
 * Return a random #gdouble equaly distributed over the range [@min..@max).
 *
 * Return value: A random number.
 **/
373 374 375 376 377 378 379 380 381 382 383 384 385
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;
}

386 387 388 389 390 391 392
/**
 * g_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 te g_random_* functions, to @seed.
 **/
393 394 395 396 397 398 399 400 401 402 403
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);
}