/**********************************************************************  random_util.c **********************************************************************  random_util - Pseudo-random number generation routines.  Copyright ?2000-2004, Stewart Adcock <stewart@linux-domain.com>  All rights reserved.  The latest version of this program should be available at:  http://gaul.sourceforge.net/  This program is free software; you can redistribute it and/or modify  it under the terms of the GNU General Public License as published by  the Free Software Foundation; either version 2 of the License, or  (at your option) any later version.  Alternatively, if your project  is incompatible with the GPL, I will probably agree to requests  for permission to use the terms of any other license.  This program is distributed in the hope that it will be useful, but  WITHOUT ANY WARRANTY WHATSOEVER.  A full copy of the GNU General Public License should be in the file  "COPYING" provided with this distribution; if not, see:  http://www.gnu.org/ **********************************************************************  Synopsis:	Portable routines for generating pseudo random numbers.		Why use these instead of the system routines?		(a) Useful interface functions. (can specify ranges or		specific distributions)		(b) System independant. (Generated numbers are		reproducible across system types.)		(c) Enables saving and restoring state.		SLang intrinsic function wrappers are provided if the		HAVE_SLANG constant is set to 1.		The algorithm I selected is the Mitchell and Moore		variant of the standard additive number generator.		The required array of numbers is populated by a		using single seed value by using a linear		congruential pseudo random number generator.		This routines have been tested and provide the same		output (within the limits of computational precision)		on (at least) the following platforms:		o Intel x86 (Linux, OpenBSD, FreeBSD)		o AMD x86-64 (Linux, FreeBSD)		o IBM Power3 (AIX)		o IBM PowerPC (Linux)		o MIPS R4K, R10K, R12K (IRIX)		o Alpha EV56 (Tru64), EV7 (Linux)		o SPARC Ultra-4 (Solaris)		This code should be thread safe.		For OpenMP code, USE_OPENMP must be defined and 		random_init() or random_seed() must be called BY A		SINGLE THREAD prior to any other function.		  References:	Standard additive number generator and the linear		congruential algorithm:		Knuth D.E., "The Art of Computer Programming",		Vol 2, 2nd ed. (1998)		General information:		Press W., Flannery B.P., Teukolsky S.A., Vetterling W.T.,		"Numerical Recipes in C: The Art of Scientific Computing"		Cambridge University Press, 2nd ed. (1992)  Usage:	o First call random_init().		o Call random_seed() to seed the PRNG (like srand()).		  Alternatively, random_tseed() to use system clock for seed.		o random_rand() is a rand() replacement, which is available		  for use but I would recommend the wrapper functions:		  random_int(), random_int_range() etc.		o random_get_state() and random_set_state() may be used		  to set, save, restore, and query the current state.		These functions can be tested by compiling with		something like:		gcc -o testrand random_util.c -DRANDOM_UTIL_TEST  To do:	Gaussian/Boltzmann distributions.		Need a proper test of randomness.		Properly implement stack of states.		I could do with some handy documentation. **********************************************************************/#include "gaul/random_util.h"/* * PRNG constants.  Don't mess with these values willie-nillie. */#define RANDOM_MM_ALPHA	55#define RANDOM_MM_BETA	24#define RANDOM_LC_ALPHA	3#define RANDOM_LC_BETA	257#define RANDOM_LC_GAMMA	RANDOM_RAND_MAX/* * Global state variable stack. * (Implemented using singly-linked list.) *//*static GSList	*state_list=NULL;*/static random_state	current_state;static boolean		is_initialised=FALSE;THREAD_LOCK_DEFINE_STATIC(random_state_lock);/********************************************************************** random_rand() Synopsis:	Replacement for the standard rand().		Returns a new pseudo-random value from the sequence, in		the range 0 to RANDOM_RAND_MAX inclusive, and updates		global state for next call.  size should be non-zero,		and state should be initialized.  parameters:  return:	none  last updated:	30/12/00 **********************************************************************/unsigned int random_rand(void)  {  unsigned int val;  if (!is_initialised) die("Neither random_init() or random_seed() have been called.");  THREAD_LOCK(random_state_lock);  val = (current_state.v[current_state.j]+current_state.v[current_state.k])        & RANDOM_RAND_MAX;  current_state.x = (current_state.x+1) % RANDOM_NUM_STATE_VALS;  current_state.j = (current_state.j+1) % RANDOM_NUM_STATE_VALS;  current_state.k = (current_state.k+1) % RANDOM_NUM_STATE_VALS;  current_state.v[current_state.x] = val;  THREAD_UNLOCK(random_state_lock);  return val;  } /**********************************************************************  random_seed()  synopsis:	Set seed for pseudo random number generator.		Uses the linear congruential algorithm to fill		state array.  parameters:	const unsigned int seed		Seed value.  return:	none  last updated: 04 May 2004 **********************************************************************/void random_seed(const unsigned int seed)  {   int	i; #if USE_OPENMP == 1  if (is_initialised == FALSE)    {    omp_init_lock(&random_state_lock);    is_initialised = TRUE;    }#else  is_initialised = TRUE;#endif  THREAD_LOCK(random_state_lock);  current_state.v[0]=(seed & RANDOM_RAND_MAX);  for(i=1; i<RANDOM_NUM_STATE_VALS; i++)    current_state.v[i] = (RANDOM_LC_ALPHA * current_state.v[i-1]                          + RANDOM_LC_BETA) & RANDOM_RAND_MAX;  current_state.j = 0;  current_state.k = RANDOM_MM_ALPHA-RANDOM_MM_BETA;  current_state.x = RANDOM_MM_ALPHA-0;  THREAD_UNLOCK(random_state_lock);  return;  } /**********************************************************************  random_tseed()  synopsis:	Set seed for pseudo random number generator from		the system time.  parameters:  return:	none  last updated:	28/01/01 **********************************************************************/void random_tseed(void)  {   random_seed((unsigned int) (time(NULL)%RANDOM_RAND_MAX));  return;  } /*************************************************************************  random_get_state()  synopsis:	Retrieve current state.  This can be used for saving		the current state.  parameters:	none  return:	random_state state  last updated:	16/05/00*************************************************************************/random_state random_get_state(void)  {  return current_state;  }/*************************************************************************  random_set_state()  synopsis:	Replace current state with specified state.		This can be used for restoring a saved state.  parameters:	random_state state  return:	none  last updated:	16/05/00*************************************************************************/void random_set_state(random_state state)  {  current_state = state;  return;  }/*************************************************************************  random_get_state_str()  synopsis:	Retrieve current state encoded as a static string.		This can be used for saving the current state.  parameters:	none  return:	char *  last updated:	28/01/01*************************************************************************/char *random_get_state_str(void)  {  return (char *) &current_state;  }/*************************************************************************  random_get_state_str_len()  synopsis:	Retrieve the length of the string encoded current state.		This can be used for saving the current state.  parameters:	none  return:	char *  last updated:	28/01/01*************************************************************************/unsigned int random_get_state_str_len(void)  {  return (unsigned int) sizeof(random_state);  }/*************************************************************************  random_set_state_str()  synopsis:	Replace current state with specified state.		This can be used for restoring a saved state.  parameters:	char *  return:	none  last updated:	28/01/01*************************************************************************/void random_set_state_str(char *state)  {  /* This causes potential unaligned trap on Alpha CPUs. */  current_state = *((random_state *)state);  return;  }/**********************************************************************  random_init()  synopsis:	Initialise random number generators.  Should be		called prior to any of these functions being used.		Seeds the PRNG with a seed of 1.  parameters:	none  return:	none  last updated:	30 May 2002 **********************************************************************/void random_init(void)  {  random_seed(1);#if RANDOM_DEBUG>0  printf("DEBUG: Random number routines initialised.\n");#endif  return;  }/**********************************************************************  random_isinit()  synopsis:	Whether these routines have been initialised.  parameters:  return:	none  last updated:	30/12/00 **********************************************************************/boolean random_isinit(void)  {  return is_initialised;  }/**********************************************************************  PRNG interface routines. **********************************************************************//**********************************************************************  random_boolean()  synopsis:	Return TRUE 50% of the time.  parameters:  return:	none  last updated:	16/05/00 **********************************************************************/boolean random_boolean(void)  {  return (random_rand() <= RANDOM_RAND_MAX/2);  }/**********************************************************************  random_boolean_prob()  synopsis:	Return TRUE (prob*100)% of the time.  parameters:  return:	TRUE or FALSE.  last updated:	16/05/00 **********************************************************************/boolean random_boolean_prob(const double prob)  {  return (random_rand() <= (unsigned int)(prob*(double)RANDOM_RAND_MAX));  }/**********************************************************************  random_int()  synopsis:	Return a random integer between 0 and (N-1) inclusive.  parameters:  return:  last updated:	05 Jun 2003 **********************************************************************/unsigned int random_int(const unsigned int max)  {/*  return (int)((double)random_rand()*max/RANDOM_RAND_MAX);*/  if (max==0)    return 0;  else    return random_rand()%max;  }/**********************************************************************  random_int_range()  synopsis:	Return a random integer between min and max-1 inclusive.  parameters:  return:  last updated:	05 Jun 2003 **********************************************************************/int random_int_range(const int min, const int max)  {/*  return (random_rand()*(max-min)/RANDOM_RAND_MAX + min);*/  if (max-min==0)    return max;  else    return min + (random_rand()%(max-min));  }