00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00066 #ifndef __random_h
00067 #define __random_h
00068 #include <math.h>
00069 #include <cstdlib>
00070 #include <complex.h>
00071
00072 #include <binary.h>
00073 namespace SPUC {
00074 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00075
00076 typedef complex<double> double_complex;
00087 class Random_Generator {
00088 public:
00090 Random_Generator();
00092 void randomize();
00094 void reset();
00096 void reset(unsigned int seed);
00098 unsigned int random_int()
00099 {
00100 if( left == 0 ) reload();
00101 --left;
00102
00103 register unsigned int s1;
00104 s1 = *pNext++;
00105 s1 ^= (s1 >> 11);
00106 s1 ^= (s1 << 7) & 0x9d2c5680U;
00107 s1 ^= (s1 << 15) & 0xefc60000U;
00108 return ( s1 ^ (s1 >> 18) );
00109 }
00111 double random_01() { return ( (double(random_int()) + 0.5) * 2.3283064365386963e-10 ); }
00113 double random_01_closed() { return ( double(random_int()) * 2.3283064370807974e-10 ); }
00114 protected:
00115 private:
00117 unsigned long lastSeed;
00119 static const unsigned int MMAXINT;
00121 static const unsigned int MAGIC;
00122
00124 static unsigned int state[624];
00126 static unsigned int *pNext;
00128 static int left;
00130 void set_seed( unsigned int oneSeed );
00132 void set_seed( unsigned int *const bigSeed );
00134 void reload()
00135 {
00136
00137
00138 register unsigned int *p = state;
00139 register int i;
00140 for( i = 624 - 397; i--; ++p )
00141 *p = twist( p[397], p[0], p[1] );
00142 for( i = 397; --i; ++p )
00143 *p = twist( p[397-624], p[0], p[1] );
00144 *p = twist( p[397-624], p[0], state[0] );
00145
00146 left = 624, pNext = state;
00147 }
00149 unsigned int hiBit( const unsigned int& u ) const { return u & 0x80000000U; }
00151 unsigned int loBit( const unsigned int& u ) const { return u & 0x00000001U; }
00153 unsigned int loBits( const unsigned int& u ) const { return u & 0x7fffffffU; }
00155 unsigned int mixBits( const unsigned int& u, const unsigned int& v ) const
00156 { return hiBit(u) | loBits(v); }
00158 unsigned int twist( const unsigned int& m, const unsigned int& s0, const unsigned int& s1 ) const
00159 { return m ^ (mixBits(s0,s1)>>1) ^ (loBit(s1) ? MAGIC : 0U); }
00160
00161 };
00162
00165
00167 void RNG_reset(unsigned long seed);
00169 void RNG_reset();
00171 void RNG_randomize();
00173
00178 class Bernoulli_RNG {
00179 public:
00181 Bernoulli_RNG(double prob) { setup(prob); }
00183 Bernoulli_RNG() { p=0.5; }
00185 void setup(double prob)
00186 {
00187
00188 p = prob;
00189 }
00191 double get_setup() const { return p; }
00193 bin operator()() { return sample(); }
00195 bin sample() { return bin( RNG.random_01() < p ? 1 : 0 ); }
00196 protected:
00197 private:
00199 double p;
00201 Random_Generator RNG;
00202 };
00203
00219 class I_Uniform_RNG {
00220 public:
00222 I_Uniform_RNG(int min=0, int max=1);
00224 void setup(int min, int max);
00226 void get_setup(int &min, int &max) const;
00228 int operator()() { return sample(); }
00230 int sample() { return ( (int)floor(RNG.random_01() * (hi - lo + 1)) + lo ); }
00231 protected:
00232 private:
00234 int lo;
00236 int hi;
00238 Random_Generator RNG;
00239 };
00240
00245 class Uniform_RNG {
00246 public:
00248 Uniform_RNG(double min=0, double max=1.0);
00250 void setup(double min, double max);
00252 void get_setup(double &min, double &max) const;
00254 double operator()() { return ( sample()* (hi_bound - lo_bound) + lo_bound ); }
00256 double sample() { return RNG.random_01(); }
00257 protected:
00258 private:
00260 double lo_bound, hi_bound;
00262 Random_Generator RNG;
00263 };
00264
00269 class Exponential_RNG {
00270 public:
00272 Exponential_RNG(double lambda=1.0);
00274 void setup(double lambda) { l=lambda; }
00276 double get_setup() const;
00278 double operator()() { return sample(); }
00279 protected:
00280 private:
00282 double sample() { return ( -log(RNG.random_01()) / l ); }
00284 double l;
00286 Random_Generator RNG;
00287 };
00288
00293 class Normal_RNG {
00294 public:
00296 Normal_RNG(double meanval, double variance) { setup(meanval, variance); };
00298 Normal_RNG() { mean = 0.0; sigma = 1.0; odd = true; }
00300 void setup(double meanval, double variance)
00301 { mean = meanval; sigma = sqrt(variance); odd = true; }
00303 void get_setup(double &meanval, double &variance) const;
00305 double operator()() { return (sigma*sample()+mean); }
00307 double sample()
00308 {
00309 double a,b;
00310
00311 if (odd) {
00312 a = TWOPI * RNG.random_01();
00313 b = sqrt(-2.0 * log(RNG.random_01()));
00314 s1 = b * cos(a);
00315 s2 = b * sin(a);
00316 odd = !odd;
00317 return s1;
00318 } else {
00319 odd = !odd;
00320 return s2;
00321 }
00322 }
00323 protected:
00324 private:
00326 double mean, sigma, s1, s2;
00328 bool odd;
00330 Random_Generator RNG;
00331 };
00332
00337 class Laplace_RNG {
00338 public:
00340 Laplace_RNG(double meanval=0.0, double variance=1.0);
00342 void setup(double meanval, double variance);
00344 void get_setup(double &meanval, double &variance) const;
00346 double operator()() { return sample(); }
00348 double sample()
00349 {
00350 u=RNG.random_01();
00351 if(u<0.5)
00352 l=log(2.0*u);
00353 else
00354 l=-log(2.0*(1-u));
00355
00356 l *= sqrt(var/2.0);
00357 l += mean;
00358
00359 return l;
00360 }
00361 protected:
00362 private:
00364 double mean, var, u, l;
00366 Random_Generator RNG;
00367 };
00368
00373 class Complex_Normal_RNG {
00374 public:
00376 Complex_Normal_RNG(double_complex mean, double variance) { setup(mean, variance); }
00378 Complex_Normal_RNG() { m = 0.0; sigma=1.0; }
00380 void setup(double_complex mean, double variance) { m = mean; sigma = sqrt(variance); }
00382 void get_setup(double_complex &mean, double &variance) { mean = m; variance = sigma*sigma; }
00384 double_complex operator()() { return sigma*sample()+m; }
00386 double_complex sample()
00387 {
00388
00389 double a = TWOPI*RNG.random_01(), b = sqrt(-log(RNG.random_01()));
00390 return double_complex(b * cos(a), b * sin(a));
00391 }
00392 protected:
00393 private:
00395 double_complex m;
00397 double sigma;
00399 Random_Generator RNG;
00400 };
00401
00406 class AR1_Normal_RNG {
00407 public:
00409 AR1_Normal_RNG(double meanval=0.0, double variance=1.0, double rho=0.0);
00411 void setup(double meanval, double variance, double rho);
00413 void get_setup(double &meanval, double &variance, double &rho) const;
00415 void reset();
00417 double operator()() { return sample(); }
00418
00419 protected:
00420 private:
00422 double sample();
00424 double my_mean, mem, r, factr, mean, var, r1, r2;
00426 bool odd;
00428 Random_Generator RNG;
00429 };
00430
00435 typedef Normal_RNG Gauss_RNG;
00436
00441 typedef AR1_Normal_RNG AR1_Gauss_RNG;
00442
00447 class Weibull_RNG {
00448 public:
00450 Weibull_RNG(double lambda=1.0, double beta=1.0);
00451
00453 void setup(double lambda, double beta);
00455 void get_setup(double &lambda, double &beta) { lambda=l; beta=b; }
00457 double operator()() { return sample(); }
00458 protected:
00459 private:
00461 double sample();
00463 double l, b;
00465 double mean, var;
00467 Random_Generator RNG;
00468 };
00469
00474 class Rayleigh_RNG {
00475 public:
00477 Rayleigh_RNG(double sigma=1.0);
00478
00480 void setup(double sigma) { sig=sigma; }
00482 double get_setup() { return sig; }
00484 double operator()() { return sample(); }
00485 protected:
00486 private:
00488 double sample();
00490 double sig;
00492 Random_Generator RNG;
00493 };
00494
00499 class Rice_RNG {
00500 public:
00502 Rice_RNG(double sigma=1.0, double _s=1.0);
00504 void setup(double sigma, double _s) { sig=sigma; s=_s; }
00506 void get_setup(double &sigma, double &_s) { sigma=sig; _s=s; }
00508 double operator()() { return sample(); }
00509 protected:
00510 private:
00512 double sample();
00514 double sig, s;
00516 Random_Generator RNG;
00517 };
00518
00521
00522
00524
00525
00526
00527 inline void Random_Generator::reset()
00528 {
00529 set_seed(lastSeed);
00530 }
00531
00532 inline void Random_Generator::reset(unsigned int seed)
00533 {
00534 lastSeed = seed;
00535 set_seed(lastSeed);
00536 }
00537
00538 inline void Random_Generator::set_seed( unsigned int oneSeed )
00539 {
00540
00541 register unsigned int *s;
00542 register int i;
00543 for( i = 624, s = state;
00544 i--;
00545 *s = oneSeed & 0xffff0000,
00546 *s++ |= ( (oneSeed *= 69069U)++ & 0xffff0000 ) >> 16,
00547 (oneSeed *= 69069U)++ ) {}
00548 reload();
00549 }
00550
00551 inline void Random_Generator::set_seed( unsigned int *const bigSeed )
00552 {
00553
00554
00555
00556
00557
00558 register unsigned int *s = state, *b = bigSeed;
00559 register int i = 624;
00560 for( ; i--; *s++ = *b++ & 0xffffffff ) {}
00561 reload();
00562 }
00563
00564 inline double AR1_Normal_RNG::sample()
00565 {
00566 double s;
00567
00568 if (odd) {
00569 r1 = RNG.random_01();
00570 r2 = RNG.random_01();
00571 s = sqrt(factr * log(r2)) * cos(TWOPI * r1);
00572 } else
00573 s = sqrt(factr * log(r2)) * sin(TWOPI * r1);
00574
00575 odd = !odd;
00576
00577 mem = s + r * mem;
00578 s = mem + mean;
00579
00580 return s;
00581 }
00582
00583 inline double Weibull_RNG::sample()
00584 {
00585 return ( pow(-log(RNG.random_01()), 1.0/b) / l );
00586 }
00587
00588 inline double Rayleigh_RNG::sample()
00589 {
00590 double r1, r2;
00591 double s1, s2, samp;
00592 double _hypot(double x, double y);
00593
00594 r1 = RNG.random_01();
00595 r2 = RNG.random_01();
00596 s1 = sqrt(-2.0 * log(r2)) * cos(TWOPI * r1);
00597 s2 = sqrt(-2.0 * log(r2)) * sin(TWOPI * r1);
00598
00599
00600 samp = sig * _hypot(s1, s2);
00601
00602 return samp;
00603 }
00604
00605 inline double Rice_RNG::sample()
00606 {
00607 double r1, r2;
00608 double s1, s2, samp;
00609 double m1 = 0.0;
00610 double m2 = sqrt(s*s - m1*m1);
00611 double _hypot(double x, double y);
00612
00613 r1 = RNG.random_01();
00614 r2 = RNG.random_01();
00615 s1 = sqrt(-2.0 * log(r2)) * cos(TWOPI * r1);
00616 s2 = sqrt(-2.0 * log(r2)) * sin(TWOPI * r1);
00617
00618
00619 samp = _hypot(sig*s1+m1, sig*s2+m2);
00620
00621 return samp;
00622 }
00623
00624
00625 bin randb(void);
00626 double randu(void);
00627 int randi(int low, int high);
00628 double randn(void);
00629 double_complex randn_c(void);
00630
00631 void RNG_reset(unsigned long seed);
00632 void RNG_reset(void);
00633 void RNG_randomize(void);
00634 #endif
00635 }
00636 #endif // __random_h