twister.h
Go to the documentation of this file.
00001 /* $Id$ */
00002 // Author: John Wu <John.Wu at ACM.org>
00003 // Copyright 2000-2012 the Regents of the University of California
00004 #ifndef IBIS_TWISTER_H
00005 #define IBIS_TWISTER_H
00006 
00019 #include <time.h>       // time_t time clock_t clock
00020 #include <limits.h>     // DBL_MIN, ULONG_MAX
00021 #include <float.h>      // DBL_MIN
00022 #include <math.h>       // sqrt, log, exp, pow
00023 
00024 #include <vector>       // std::vector<double> used by discrateZip1
00025 
00026 namespace ibis {
00027     class uniformRandomNumber;  // the abstract base class
00028     class MersenneTwister;      // concrete uniform random number generator
00029     class discretePoisson;      // discrete Poisson
00030     class discretePoisson1;     // the spcial case for exp(-x)
00031     class discreteZipf;         // discrete Zipf 1/x^a (a > 0)
00032     class discreteZipf1;        // 1/x
00033     class discreteZipf2;        // 1/x^2
00034     class randomGaussian;       // continuous Gaussian
00035     class randomPoisson;        // continuous Posson
00036     class randomZipf;           // continuous Zipf
00037 };
00038 
00040 class ibis::uniformRandomNumber {
00041 public:
00042     virtual double operator()() = 0;
00043 };
00044 
00047 class ibis::MersenneTwister : public ibis::uniformRandomNumber {
00048 public:
00051     MersenneTwister() {setSeed(time(0) ^ clock());}
00053     MersenneTwister(unsigned seed) {setSeed(seed);}
00054 
00056     virtual double operator()() {return nextDouble();}
00058     int nextInt() {return next();}
00059     long nextLong() {return next();}
00060     float nextFloat() {return 2.3283064365386962890625e-10*next();}
00061     double nextDouble() {return 2.3283064365386962890625e-10*next();}
00063     unsigned next(unsigned r) {return static_cast<unsigned>(r*nextDouble());}
00064 
00066     void setSeed(unsigned seed) {
00067         for (int i=0; i<624; i++) {
00068             mt[i] = seed & 0xffff0000;
00069             seed = 69069 * seed + 1;
00070             mt[i] |= (seed & 0xffff0000) >> 16;
00071             seed = 69069 * seed + 1;
00072         }
00073         mti = 624;
00074     }
00075 
00077     unsigned next() {
00078         unsigned y;
00079         if (mti >= 624) { /* generate 624 words at one time */
00080             static unsigned mag01[2]={0x0, 0x9908b0df};
00081             int kk;
00082 
00083             for (kk=0; kk<227; kk++) {
00084                 y = (mt[kk]&0x80000000)|(mt[kk+1]&0x7fffffff);
00085                 mt[kk] = mt[kk+397] ^ (y >> 1) ^ mag01[y & 0x1];
00086             }
00087             for (; kk<623; kk++) {
00088                 y = (mt[kk]&0x80000000)|(mt[kk+1]&0x7fffffff);
00089                 mt[kk] = mt[kk-227] ^ (y >> 1) ^ mag01[y & 0x1];
00090             }
00091             y = (mt[623]&0x80000000)|(mt[0]&0x7fffffff);
00092             mt[623] = mt[396] ^ (y >> 1) ^ mag01[y & 0x1];
00093             mti = 0;
00094         }
00095   
00096         y = mt[mti++];
00097         y ^= (y >> 11);
00098         y ^= (y << 7) & 0x9d2c5680;
00099         y ^= (y << 15) & 0xefc60000;
00100         y ^= (y >> 18);
00101 
00102         return y; 
00103     }
00104 
00105 private:
00106     // all data members are private
00107     int mti;
00108     unsigned mt[624]; /* the array for the state vector  */
00109 }; // class MersenneTwister
00110 
00111 // ********** the following random number generators need a uniformation
00112 // ********** random number generator as the input
00113 
00115 class ibis::randomPoisson {
00116 public:
00118     randomPoisson(uniformRandomNumber& ur) : urand(ur) {}
00119     double operator()() {return next();}
00120     double next() {return -log(urand());}
00121 
00122 private:
00123     uniformRandomNumber& urand;
00124 }; // ibis::randomPoisson
00125 
00129 class ibis::randomGaussian {
00130 public:
00132     randomGaussian(uniformRandomNumber& ur)
00133         : urand(ur), has_extra(false), extra(0.0) {}
00135     double operator()() {return next();}
00137     double next() {
00138         if (has_extra) { /* has extra value from the previous run */
00139             has_extra = false;
00140             return extra;
00141         }
00142         else { /* Box-Mueller transformation */
00143             double v1, v2, r, fac;
00144             do {
00145                 v1 = 2.0 * urand() - 1.0;
00146                 v2 = 2.0 * urand() - 1.0;
00147                 r = v1 * v1 + v2 * v2;
00148             } while (r >= 1.0 || r <= 0.0);
00149             fac = sqrt((-2.0 * log(r))/r);
00150             has_extra = false;
00151             extra = v2 * fac;
00152             v1 *= fac;
00153             return v1;
00154         }
00155     }
00156 
00157 private:
00158     uniformRandomNumber& urand;
00159     bool has_extra;
00160     double extra;
00161 }; // ibis::randomGaussian
00162 
00165 class ibis::randomZipf {
00166 public:
00168     randomZipf(uniformRandomNumber& ur, double a=1) : urand(ur), alpha(a-1) {}
00170     double operator()() {return next();}
00172     double next() {
00173         if (alpha > 0.0)
00174             return (exp(-log(1 - urand())/(alpha)) - 1);
00175         else
00176             return (1.0 / (1.0 - urand()) - 1.0);
00177     }
00178 
00179 private:
00180     uniformRandomNumber& urand;
00181     const double alpha; // Zipf exponent - 1
00182 }; // ibis::randomZipf
00183 
00186 class ibis::discretePoisson {
00187 public:
00188     discretePoisson(ibis::uniformRandomNumber& ur,
00189                     const double lam=1.0, long m=0)
00190         : min0(m), lambda(lam), urand(ur) {init();}
00191 
00192     long operator()() {return next();}
00193     long next() {
00194         long k;
00195         double u, x;
00196         while (true) {
00197             u = ym * (urand)();
00198             x = - lambda * log(u * laminv);
00199             k = static_cast<long>(x + 0.5);
00200             if (k <= k0 && k-x <= xm)
00201                 return k;
00202             else if (u >= -exp(laminv*k+laminv2)*lambda-exp(laminv*k))
00203                 return k;
00204         }
00205     } // next integer random number
00206 
00207 private:
00208     // private member variables
00209     long min0, k0;
00210     double lambda, laminv, laminv2, xm, ym;
00211     uniformRandomNumber& urand;
00212 
00213     // private functions
00214     void init() { // check input parameters and initialize three constants
00215         if (! (lambda > DBL_MIN))
00216             lambda = 1.0;
00217         laminv = -1.0 / lambda;
00218         laminv2 = 0.5*laminv;
00219         k0 = static_cast<long>(1.0+min0+1.0/(1.0-exp(laminv)));
00220         ym = -exp((min0+0.5)*laminv)*lambda - exp(min0*laminv);
00221         xm = min0 - log(ym*laminv);
00222     }
00223 }; // class discretePoisson
00224 
00226 class ibis::discretePoisson1 {
00227 public:
00228     discretePoisson1(ibis::uniformRandomNumber& ur) : urand(ur) {init();}
00229 
00230     long operator()() {return next();}
00231     long next() {
00232         long k;
00233         double u, x;
00234         while (true) {
00235             u = ym * (urand)();
00236             x = - log(-u);
00237             k = static_cast<long>(x + 0.5);
00238             if (k <= k0 && k-x <= xm)
00239                 return k;
00240             else if (u >= -exp(-static_cast<double>(k)-0.5) -
00241                      exp(-static_cast<double>(k)))
00242                 return k;
00243         }
00244     } // next integer random number
00245 
00246 private:
00247     // private member variables
00248     double xm, ym;
00249     long k0;
00250     uniformRandomNumber& urand;
00251 
00252     // private functions
00253     void init() { // check input parameters and initialize three constants
00254         k0 = static_cast<long>(1.0+1.0/(1.0-exp(-1.0)));
00255         ym = - exp(-0.5) - 1.0;
00256         xm = - log(-ym);
00257     }
00258 }; // class discretePoisson1
00259 
00264 class ibis::discreteZipf {
00265 public:
00266     discreteZipf(ibis::uniformRandomNumber& ur, double a=2.0,
00267                  unsigned long imax = ULONG_MAX)
00268         : urand(ur), max0(imax), alpha(a) {init();}
00269 
00271     unsigned long operator()() {return next();}
00272     unsigned long next() {
00273         if (alpha > 1.0) { // rejection-inversion
00274             while (true) {
00275                 double ur = (urand)();
00276                 ur = hxm + ur * hx0;
00277                 double x = Hinv(ur);
00278                 unsigned long k = static_cast<unsigned long>(0.5+x);
00279                 if (k - x <= ss)
00280                     return k;
00281                 else if (ur >= H(0.5+k) - exp(-log(k+1.0)*alpha))
00282                     return k;
00283             }
00284         }
00285         else { // simple rejection
00286             unsigned long k = ((long) (urand() * max0)) % max0;
00287             double freq = pow((1.0+k), -alpha);
00288             while (urand() >= freq) {
00289                 k = ((long) (urand() * max0)) % max0;
00290                 freq = pow((1.0+k), -alpha);
00291             }
00292             return k;
00293         }
00294     } // next
00295 
00296 private:
00297     // private member variables
00298     uniformRandomNumber& urand;
00299     long unsigned max0;
00300     double alpha, alpha1, alphainv, hx0, hxm, ss;
00301 
00302     // private member function
00303     double H(double x) {return (exp(alpha1*log(1.0+x)) * alphainv);}
00304     double Hinv(double x) {return exp(alphainv*log(alpha1*x)) - 1.0;}
00305     void init() {
00306         // enforce the condition that alpha >= 0 and max0 > 1
00307         if (max0 <= 1)
00308             max0 = 100;
00309         if (! (alpha >= 0.0))
00310             alpha = 1.0;
00311         if (alpha > 1.0) {
00312             // the rejection-inversion algorithm of W. Hormann and
00313             // G. Derflinger
00314             alpha1 = 1.0 - alpha;
00315             alphainv = 1.0 / alpha1;
00316             hxm = H(max0 + 0.5);
00317             hx0 = H(0.5) - 1.0 - hxm;
00318             ss = 1 - Hinv(H(1.5)-exp(-alpha*log(2.0)));
00319         }
00320         else { // use a simple rejection scheme
00321             alpha1 = 0.0;
00322             alphainv = 0.0;
00323             hxm = 0.0;
00324             hx0 = 0.0;
00325             ss  = 0.0;
00326         }
00327     }
00328 }; // Zipf distribution
00329 
00332 class ibis::discreteZipf2 {
00333 public:
00334     discreteZipf2(ibis::uniformRandomNumber& ur,
00335                   unsigned long imax = ULONG_MAX) :
00336         max0(imax), urand(ur) {init();}
00337 
00339     unsigned long operator()() {return next();}
00340     unsigned long next() {
00341         while (true) {
00342             double ur = (urand)();
00343             ur = hxm + ur * hx0;
00344             double x = Hinv(ur);
00345             unsigned long k = static_cast<unsigned long>(0.5+x);
00346             if (k - x <= ss)
00347                 return k;
00348             else if (ur >= H(0.5+k) - 1.0/((1.0+x)*(1.0+x)))
00349                 return k;
00350         }
00351     } // next
00352 
00353 private:
00354     // private member variables
00355     double hx0, hxm, ss;
00356     long unsigned max0;
00357     uniformRandomNumber& urand;
00358 
00359     // private member function
00360     double H(double x) {return -1.0 / (1.0 + x);}
00361     double Hinv(double x) {return (- 1.0 / x) - 1.0;}
00362     void init() {
00363         hxm = H(max0+0.5);
00364         hx0 = - 5.0/3.0 - hxm;
00365         ss = 1 - Hinv(H(1.5)-0.25);
00366     }
00367 }; // Zipf2 distribution
00368 
00372 class ibis::discreteZipf1 {
00373 public:
00374     discreteZipf1(ibis::uniformRandomNumber& ur, unsigned long imax = 100) :
00375         card(imax+1), cpd(imax+1), urand(ur) {init();}
00376 
00378     unsigned long operator()() {return next();}
00379     unsigned long next() {
00380         double ur = (urand)();
00381         if (ur <= cpd[0]) return 0;
00382         // return the minimal i such that cdf[i] >= ur
00383         unsigned long i, j, k;
00384         i = 0;
00385         j = card-1;
00386         k = (i + j) / 2;
00387         while (i < k) {
00388             if (cpd[k] > ur)
00389                 j = k;
00390             else if (cpd[k] < ur)
00391                 i = k;
00392             else
00393                 return k;
00394             k = (i + j) / 2;
00395         }
00396         if (cpd[i] >= ur)
00397             return i;
00398         else
00399             return j;
00400     } // next
00401 
00402 private:
00403     // private member variables
00404     const unsigned long card;
00405     std::vector<double> cpd; // cumulative probability distribution
00406     uniformRandomNumber& urand;
00407 
00408     // private member function
00409     void init() { // generates the cpd
00410         const unsigned n = cpd.size();
00411         if (n < 2 || n > 1024*1024 || card != n)
00412             throw "imax must be in [2, 1000000]";
00413 
00414         cpd[0] = 1.0;
00415         for (unsigned i = 1; i < n; ++i)
00416             cpd[i] = cpd[i-1] + 1.0 / (1.0 + i);
00417         double ss = 1.0 / cpd.back();
00418         for (unsigned i = 0; i < n; ++i)
00419             cpd[i] *= ss;
00420     } // init
00421 }; // Zipf1 distribution
00422 #endif // IBIS_TWISTER_H

Make It A Bit Faster
Contact us
Disclaimers
FastBit source code
FastBit mailing list archive