📄 garandom.c
字号:
// $Header: /usr/people/mbwall/src/galib/ga/RCS/garandom.C,v 1.8 1999/03/27 19:25:44 mbwall Exp $
/* ----------------------------------------------------------------------------
random.C
mbwall 5sep95
Copyright (c) 1995 Massachusetts Institute of Technology
all rights reserved
DESCRIPTION:
Random number stuff for use in GAlib.
---------------------------------------------------------------------------- */
#include <ga/garandom.h>
#include <time.h>
#include <math.h>
#include <string.h>
static void bitseed(unsigned int seed=1);
// If the machine has multiple processes, use the PID to help make the random
// number generator seed more random.
#if defined(USE_PID)
#include <unistd.h>
#define _GA_PID * getpid()
#else
#define _GA_PID
#endif
// Return a string indicating which random number generator was compiled-in to
// the library.
const char*
GAGetRNG() {
#if defined(USE_RAN1)
return "RAN1";
#elif defined(USE_RAN2)
return "RAN2";
#elif defined(USE_RAN3)
return "RAN3";
#elif defined(USE_RAND)
return "RAND";
#elif defined(USE_RANDOM)
return "RANDOM";
#elif defined(USE_RAND48)
return "RAND48";
#else
return "UNKNOWN";
#endif
}
// Seed the random number generator with an appropriate value. We seed both
// the random number generator and the random bit generator. Set the seed only
// if a seed is not specified. If a seed is specified, then set the seed to
// the specified value and use it. We remember the seed so that multiple calls
// to this function with the same seed do not reset the generator. Subsequent
// calls to this function with a different seed will initialize the generator
// to the new seed. Multiple calls with a value of 0 do nothing (we do *not*
// re-seed the generator because 0 is the default value and we don't want
// people to re-seed the generator inadvertantly).
// Some systems return a long as the return value for time, so we need to be
// sure to get whatever variation from it that we can since our seed is only an
// unsigned int.
static unsigned int seed=0;
unsigned int GAGetRandomSeed() { return seed; }
void
GARandomSeed(unsigned int s) {
if(s == 0 && seed == 0) {
unsigned long int tmp;
while(seed == 0) {
tmp = time(NULL) _GA_PID;
for(unsigned int i=0; i<BITS_IN_WORD*sizeof(unsigned int); i++)
seed += (tmp & (1 << i));
}
_GA_RND_SEED (seed);
bitseed(seed);
}
else if(s != 0 && seed != s) {
seed = s;
_GA_RND_SEED (seed);
bitseed(seed);
}
}
// Similar to setting the random seed, but this one sets it as long as the
// specified seed is non-zero.
void
GAResetRNG(unsigned int s) {
if(s != 0) {
seed = s;
_GA_RND_SEED (seed);
bitseed(seed);
}
}
// Return a number from a unit Gaussian distribution. The mean is 0 and the
// standard deviation is 1.0.
// First we generate two uniformly random variables inside the complex unit
// circle. Then we transform these into Gaussians using the Box-Muller
// transformation. This method is described in Numerical Recipes in C
// ISBN 0-521-43108-5 at http://world.std.com/~nr
// When we find a number, we also find its twin, so we cache that here so
// that every other call is a lookup rather than a calculation. (I think GNU
// does this in their implementations as well, but I don't remember for
// certain.)
double
GAUnitGaussian(){
static GABoolean cached=gaFalse;
static double cachevalue;
if(cached == gaTrue){
cached = gaFalse;
return cachevalue;
}
double rsquare, factor, var1, var2;
do{
var1 = 2.0 * GARandomDouble() - 1.0;
var2 = 2.0 * GARandomDouble() - 1.0;
rsquare = var1*var1 + var2*var2;
} while(rsquare >= 1.0 || rsquare == 0.0);
double val = -2.0 * log(rsquare) / rsquare;
if(val > 0.0) factor = sqrt(val);
else factor = 0.0; // should not happen, but might due to roundoff
cachevalue = var1 * factor;
cached = gaTrue;
return (var2 * factor);
}
// This is the random bit generator Method II from numerical recipes in C. The
// seed determines where in the cycle of numbers the generator will start, so
// we don't need full 'long' precision in the argument to the seed function.
#define IB1 1
#define IB2 2
#define IB5 16
#define IB18 131072L
#define MASK (IB1+IB2+IB5)
static unsigned long iseed;
void
bitseed(unsigned int seed) {
iseed = seed;
}
int
GARandomBit() {
if (iseed & IB18) {
iseed=((iseed ^ MASK) << 1) | IB1;
return 1;
} else {
iseed <<= 1;
return 0;
}
}
#undef MASK
#undef IB18
#undef IB5
#undef IB2
#undef IB1
// The following random number generators are from Numerical Recipes in C.
// I have split them into a seed function and random number function.
// The ran1 pseudo-random number generator. This one is OK to use as long as
// you don't call it more than about 10^8 times, so for any long GA runs you'd
// better use something with a longer period.
#if defined(USE_RAN1)
#define IA 16807L
#define IM 2147483647L
#define AM (1.0/IM)
#define IQ 127773L
#define IR 2836L
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)
static long iy=0;
static long iv[NTAB];
static long idum=0;
void
sran1(unsigned seed) {
int j;
long k;
idum = seed;
if (idum < 1) idum=1;
for (j=NTAB+7;j>=0;j--) {
k=(idum)/IQ;
idum=IA*(idum-k*IQ)-IR*k;
if (idum < 0) idum += IM;
if (j < NTAB) iv[j] = idum;
}
iy=iv[0];
}
float ran1() {
int j;
long k;
float temp;
k=(idum)/IQ;
idum=IA*(idum-k*IQ)-IR*k;
if (idum < 0) idum += IM;
j=iy/NDIV;
iy=iv[j];
iv[j] = idum;
if ((temp=AM*iy) > RNMX) return RNMX;
else return temp;
}
#undef IA
#undef IM
#undef AM
#undef IQ
#undef IR
#undef NTAB
#undef NDIV
#undef EPS
#undef RNMX
#endif
// The ran2 pseudo-random number generator. It has a period of 2 * 10^18 and
// returns a uniform random deviate on the interval (0.0, 1.0) excluding the
// end values. idum initializes the sequence, so we create a separate seeding
// function to set the seed. If you reset the seed then you re-initialize the
// sequence.
#if defined(USE_RAN2)
#define IM1 2147483563L
#define IM2 2147483399L
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014L
#define IA2 40692L
#define IQ1 53668L
#define IQ2 52774L
#define IR1 12211L
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)
static long idum2=123456789;
static long iy=0;
static long iv[NTAB];
static long idum=0;
void
sran2(unsigned int seed) {
int j;
long k;
idum = STA_CAST(long,seed);
if (idum < 1) idum=1;
idum2=(idum);
for (j=NTAB+7;j>=0;j--) {
k=(idum)/IQ1;
idum=IA1*(idum-k*IQ1)-k*IR1;
if (idum < 0) idum += IM1;
if (j < NTAB) iv[j] = idum;
}
iy=iv[0];
}
float
ran2() {
int j;
long k;
float temp;
k=(idum)/IQ1;
idum=IA1*(idum-k*IQ1)-k*IR1;
if (idum < 0) idum += IM1;
k=idum2/IQ2;
idum2=IA2*(idum2-k*IQ2)-k*IR2;
if (idum2 < 0) idum2 += IM2;
j=iy/NDIV;
iy=iv[j]-idum2;
iv[j] = idum;
if (iy < 1) iy += IMM1;
if ((temp=AM*iy) > RNMX) return RNMX;
else return temp;
}
#undef IM1
#undef IM2
#undef AM
#undef IMM1
#undef IA1
#undef IA2
#undef IQ1
#undef IQ2
#undef IR1
#undef IR2
#undef NTAB
#undef NDIV
#undef EPS
#undef RNMX
#endif
#if defined(USE_RAN3)
// The ran3 pseudo-random number generator. It is *not* linear congruential.
#define MBIG 1000000000
#define MSEED 161803398
#define MZ 0
#define FAC (1.0/MBIG)
static int inext,inextp;
static long ma[56];
void
sran3(unsigned int seed) {
long idum = seed;
long mj,mk;
int i,ii,k;
mj=MSEED-idum;
mj %= MBIG;
ma[55]=mj;
mk=1;
for (i=1;i<=54;i++) {
ii=(21*i) % 55;
ma[ii]=mk;
mk=mj-mk;
if (mk < MZ) mk += MBIG;
mj=ma[ii];
}
for (k=1;k<=4;k++)
for (i=1;i<=55;i++) {
ma[i] -= ma[1+(i+30) % 55];
if (ma[i] < MZ) ma[i] += MBIG;
}
inext=0;
inextp=31;
}
float
ran3() {
long mj;
int i,ii,k;
if (++inext == 56) inext=1;
if (++inextp == 56) inextp=1;
mj=ma[inext]-ma[inextp];
if (mj < MZ) mj += MBIG;
ma[inext]=mj;
return mj*FAC;
}
#undef MBIG
#undef MSEED
#undef MZ
#undef FAC
#endif
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -