00001
00002
00003
00004 #ifndef IBIS_TWISTER_H
00005 #define IBIS_TWISTER_H
00006
00020 #include <time.h>
00021 #include <math.h>
00022 #include <limits.h>
00023 #include <float.h>
00024
00025 #include <vector>
00026
00027 namespace ibis {
00028 class uniformRandomNumber;
00029 class MersenneTwister;
00030 class discretePoisson;
00031 class discretePoisson1;
00032 class discreteZipf;
00033 class discreteZipf1;
00034 class discreteZipf2;
00035 };
00036
00038 class ibis::uniformRandomNumber {
00039 public:
00040 virtual double operator()() = 0;
00041 };
00042
00044 class ibis::MersenneTwister : public ibis::uniformRandomNumber {
00045 public:
00046
00047
00048 MersenneTwister() {setSeed(time(0) ^ clock());}
00049 MersenneTwister(unsigned seed) {setSeed(seed);}
00050
00051 virtual double operator()() {return nextDouble();}
00052 int nextInt() {return next();}
00053 long nextLong() {return next();}
00054 float nextFloat() {return 2.3283064365386962890625e-10*next();}
00055 double nextDouble() {return 2.3283064365386962890625e-10*next();}
00056
00057 unsigned next(unsigned r) {return static_cast<unsigned>(r*nextDouble());}
00058
00059
00060
00061 double nextPoisson() {return -log(nextDouble());}
00062
00063
00064
00065 double nextZipf() {return (1.0 / (1.0 - nextDouble()) - 1.0);}
00066
00067
00068
00069 double nextZipf(double a) {
00070 if (a > 1.0)
00071 return (exp(-log(1-nextDouble())/(a-1)) - 1);
00072 else
00073 return -1.0;
00074 }
00075
00076
00077
00078 double nextGaussian() {
00079 if (has_gaussian != 0) {
00080 has_gaussian = 0;
00081 return gaussian_extra;
00082 }
00083 else {
00084 double v1, v2, r, fac;
00085 do {
00086 v1 = 4.656612873077392578125e-10*next() - 1.0;
00087 v2 = 4.656612873077392578125e-10*next() - 1.0;
00088 r = v1 * v1 + v2 * v2;
00089 } while (r >= 1.0);
00090 fac = sqrt(-2.0 * log((double) r)/r);
00091 gaussian_extra = v2 * fac;
00092 v1 *= fac;
00093 has_gaussian = 1;
00094 return v1;
00095 }
00096 }
00097
00098
00099 void setSeed(unsigned seed) {
00100 for (int i=0; i<624; i++) {
00101 mt[i] = seed & 0xffff0000;
00102 seed = 69069 * seed + 1;
00103 mt[i] |= (seed & 0xffff0000) >> 16;
00104 seed = 69069 * seed + 1;
00105 }
00106 mti = 624;
00107 has_gaussian = 0;
00108 }
00109
00110
00111 unsigned next() {
00112 unsigned y;
00113 if (mti >= 624) {
00114 static unsigned mag01[2]={0x0, 0x9908b0df};
00115 int kk;
00116
00117 for (kk=0; kk<227; kk++) {
00118 y = (mt[kk]&0x80000000)|(mt[kk+1]&0x7fffffff);
00119 mt[kk] = mt[kk+397] ^ (y >> 1) ^ mag01[y & 0x1];
00120 }
00121 for (; kk<623; kk++) {
00122 y = (mt[kk]&0x80000000)|(mt[kk+1]&0x7fffffff);
00123 mt[kk] = mt[kk-227] ^ (y >> 1) ^ mag01[y & 0x1];
00124 }
00125 y = (mt[623]&0x80000000)|(mt[0]&0x7fffffff);
00126 mt[623] = mt[396] ^ (y >> 1) ^ mag01[y & 0x1];
00127 mti = 0;
00128 }
00129
00130 y = mt[mti++];
00131 y ^= (y >> 11);
00132 y ^= (y << 7) & 0x9d2c5680;
00133 y ^= (y << 15) & 0xefc60000;
00134 y ^= (y >> 18);
00135
00136 return y;
00137 }
00138
00139 private:
00140
00141 int mti;
00142 int has_gaussian;
00143 double gaussian_extra;
00144 unsigned mt[624];
00145 };
00146
00147
00148
00149
00152 class ibis::discretePoisson {
00153 public:
00154 discretePoisson(ibis::uniformRandomNumber* ur,
00155 const double lam=1.0, long m=0)
00156 : min0(m), lambda(lam), urand(ur) {init();}
00157
00158 long operator()() {return next();}
00159 long next() {
00160 long k;
00161 double u, x;
00162 while (true) {
00163 u = ym * (*urand)();
00164 x = - lambda * log(u * laminv);
00165 k = static_cast<long>(x + 0.5);
00166 if (k <= k0 && k-x <= xm)
00167 return k;
00168 else if (u >= -exp(laminv*k+laminv2)*lambda-exp(laminv*k))
00169 return k;
00170 }
00171 }
00172
00173 private:
00174
00175 long min0, k0;
00176 double lambda, laminv, laminv2, xm, ym;
00177 uniformRandomNumber* urand;
00178
00179
00180 void init() {
00181 if (! (lambda > DBL_MIN))
00182 lambda = 1.0;
00183 laminv = -1.0 / lambda;
00184 laminv2 = 0.5*laminv;
00185 k0 = static_cast<long>(1.0+min0+1.0/(1.0-exp(laminv)));
00186 ym = -exp((min0+0.5)*laminv)*lambda - exp(min0*laminv);
00187 xm = min0 - log(ym*laminv);
00188 }
00189 };
00190
00192 class ibis::discretePoisson1 {
00193 public:
00194 discretePoisson1(ibis::uniformRandomNumber* ur) : urand(ur) {init();}
00195
00196 long operator()() {return next();}
00197 long next() {
00198 long k;
00199 double u, x;
00200 while (true) {
00201 u = ym * (*urand)();
00202 x = - log(-u);
00203 k = static_cast<long>(x + 0.5);
00204 if (k <= k0 && k-x <= xm)
00205 return k;
00206 else if (u >= -exp(-static_cast<double>(k)-0.5) -
00207 exp(-static_cast<double>(k)))
00208 return k;
00209 }
00210 }
00211
00212 private:
00213
00214 double xm, ym;
00215 long k0;
00216 uniformRandomNumber* urand;
00217
00218
00219 void init() {
00220 k0 = static_cast<long>(1.0+1.0/(1.0-exp(-1.0)));
00221 ym = - exp(-0.5) - 1.0;
00222 xm = - log(-ym);
00223 }
00224 };
00225
00230 class ibis::discreteZipf {
00231 public:
00232 discreteZipf(ibis::uniformRandomNumber* ur, double a=2.0,
00233 unsigned long v=1, unsigned long imax = ULONG_MAX) :
00234 min0(v), max0(imax), alpha(a), urand(ur) {init();}
00235
00236
00237 unsigned long operator()() {return next();}
00238 unsigned long next() {
00239 while (true) {
00240 double ur = (*urand)();
00241 ur = hxm + ur * hx0;
00242 double x = Hinv(ur);
00243 unsigned long k = static_cast<unsigned long>(0.5+x);
00244 if (k - x <= ss)
00245 return k;
00246 else if (ur >= H(0.5+k) - exp(-log(static_cast<double>
00247 (k+min0))*alpha))
00248 return k;
00249 }
00250 }
00251
00252 private:
00253
00254 long unsigned min0, max0;
00255 double alpha, alpha1, alphainv, hx0, hxm, ss;
00256 uniformRandomNumber* urand;
00257
00258
00259 double H(double x) {return (exp(alpha1*log(min0+x)) * alphainv);}
00260 double Hinv(double x) {return exp(alphainv*log(alpha1*x)) - min0;}
00261 void init() {
00262
00263 if (! (alpha > 1.0))
00264 alpha = 2.0;
00265 if (min0 < 1)
00266 min0 = 1;
00267 alpha1 = 1.0 - alpha;
00268 alphainv = 1.0 / alpha1;
00269 hxm = H(max0+0.5);
00270 hx0 = H(0.5) - exp(log(static_cast<double>(min0))*(-alpha)) - hxm;
00271 ss = 1 - Hinv(H(1.5)-exp(-alpha*log(static_cast<double>(min0)+1.0)));
00272 }
00273 };
00274
00277 class ibis::discreteZipf2 {
00278 public:
00279 discreteZipf2(ibis::uniformRandomNumber* ur,
00280 unsigned long imax = ULONG_MAX) :
00281 max0(imax), urand(ur) {init();}
00282
00284 unsigned long operator()() {return next();}
00285 unsigned long next() {
00286 while (true) {
00287 double ur = (*urand)();
00288 ur = hxm + ur * hx0;
00289 double x = Hinv(ur);
00290 unsigned long k = static_cast<unsigned long>(0.5+x);
00291 if (k - x <= ss)
00292 return k;
00293 else if (ur >= H(0.5+k) - 1.0/((1.0+x)*(1.0+x)))
00294 return k;
00295 }
00296 }
00297
00298 private:
00299
00300 double hx0, hxm, ss;
00301 long unsigned max0;
00302 uniformRandomNumber* urand;
00303
00304
00305 double H(double x) {return -1.0 / (1.0 + x);}
00306 double Hinv(double x) {return (- 1.0 / x) - 1.0;}
00307 void init() {
00308 hxm = H(max0+0.5);
00309 hx0 = - 5.0/3.0 - hxm;
00310 ss = 1 - Hinv(H(1.5)-0.25);
00311 }
00312 };
00313
00316 class ibis::discreteZipf1 {
00317 public:
00318 discreteZipf1(ibis::uniformRandomNumber* ur, unsigned long imax = 100) :
00319 card(imax+1), cpd(imax+1), urand(ur) {init();}
00320
00321
00322 unsigned long operator()() {return next();}
00323 unsigned long next() {
00324 double ur = (*urand)();
00325 if (ur <= cpd[0]) return 0;
00326
00327 unsigned long i, j, k;
00328 i = 0;
00329 j = card-1;
00330 k = (i + j) / 2;
00331 while (i < k) {
00332 if (cpd[k] > ur)
00333 j = k;
00334 else if (cpd[k] < ur)
00335 i = k;
00336 else
00337 return k;
00338 k = (i + j) / 2;
00339 }
00340 if (cpd[i] >= ur)
00341 return i;
00342 else
00343 return j;
00344 }
00345
00346 private:
00347
00348 const unsigned long card;
00349 std::vector<double> cpd;
00350 uniformRandomNumber* urand;
00351
00352
00353 void init() {
00354 const unsigned n = cpd.size();
00355 if (n < 2 || n > 1024*1024 || card != n)
00356 throw "imax must be in [2, 1000000]";
00357
00358 cpd[0] = 1.0;
00359 for (unsigned i = 1; i < n; ++i)
00360 cpd[i] = cpd[i-1] + 1.0 / (1.0 + i);
00361 double ss = 1.0 / cpd.back();
00362 for (unsigned i = 0; i < n; ++i)
00363 cpd[i] *= ss;
00364 }
00365 };
00366 #endif // IBIS_TWISTER_H