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
![]() |