Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

random.h

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*
00002  *                                   IT++                                    *
00003  *---------------------------------------------------------------------------*
00004  * Copyright (c) 1995-2002 by Tony Ottosson, Thomas Eriksson, Pål Frenger,   *
00005  * Tobias Ringström, and Jonas Samuelsson.                                   *
00006  *                                                                           *
00007  * Permission to use, copy, modify, and distribute this software and its     *
00008  * documentation under the terms of the GNU General Public License is hereby *
00009  * granted. No representations are made about the suitability of this        *
00010  * software for any purpose. It is provided "as is" without expressed or     *
00011  * implied warranty. See the GNU General Public License for more details.    *
00012  *---------------------------------------------------------------------------*/
00013 // 
00014 //  SPUC - Signal processing using C++ - A DSP library
00066 #ifndef __random_h
00067 #define __random_h
00068 #include <math.h>
00069 #include <cstdlib>
00070 #include <complex.h>
00071 //#include "spucassert.h"
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 ); } //*2^-32+2^-33
00113   double random_01_closed() { return ( double(random_int()) * 2.3283064370807974e-10 ); } // * 1/(2^32-1)
00114  protected:
00115  private:
00117   unsigned long lastSeed;
00119   static const unsigned int MMAXINT;   // largest value from randInt()
00121   static const unsigned int MAGIC;    // magic constant
00122   
00124   static unsigned int state[624];     // internal state
00126   static unsigned int *pNext;         // next value to get from state
00128   static int left;                    // number of values left before reload needed
00130   void set_seed( unsigned int oneSeed );
00132   void set_seed( unsigned int *const bigSeed );
00134   void reload()
00135     {
00136       // Generate N new values in state
00137       // Made clearer and faster by Matthew Bellew (matthew.bellew@home.com)
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   //  static unsigned int hash( time_t t, clock_t c );
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 //      spuc_assert(prob>=0.0 && prob<=1.0, "The bernoulli source probability must be between 0 and 1");
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       // Make a const of 2.0*M_PI
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 // -------------------- INLINES ----------------------------------------------
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         // Seed the generator with a simple uint32
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)++ ) {}  // hard to read, but fast
00548         reload();
00549 }
00550 
00551 inline void Random_Generator::set_seed( unsigned int *const bigSeed )
00552 {
00553         // Seed the generator with an array of 624 uint32's
00554         // There are 2^19937-1 possible initial states.  This function
00555         // allows any one of those to be chosen by providing 19937 bits.
00556         // Theoretically, the array can contain any values except all zeroes.
00557         // Just call seed() if you want to get array from /dev/urandom
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     // s1 and s2 are N(0,1) and independent
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     // s1 and s2 are N(0,1) and independent
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 // Set the seed of the Global Random Number Generator
00631 void RNG_reset(unsigned long seed);
00632 void RNG_reset(void);
00633 void RNG_randomize(void);
00634 #endif
00635 }
00636 #endif // __random_h

Generated on Fri Sep 16 11:02:28 2005 for spuc by  doxygen 1.4.4