Random.cpp
00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include "Random.h"
00043 #include "MathConstant.h"
00044 #include "Error.h"
00045
00046 #include <math.h>
00047
00048
00049 #define MERSENNE_N 624
00050 #define MERSENNE_M 397
00051
00052
00053 #define MATRIX_A 0x9908b0dfUL
00054
00055
00056 #define UPPER_MASK 0x80000000UL
00057
00058
00059 #define LOWER_MASK 0x7fffffffUL
00060
00061
00062
00063 #define IA 16807
00064 #define IM 2147483647
00065 #define AM (1.0 / IM)
00066 #define IQ 127773
00067 #define IR 2836
00068 #define NTAB 32
00069 #define NDIV (1+(IM-1)/NTAB)
00070 #define RNMX (1.0-EPS)
00071
00072 Random::Random(long s)
00073 {
00074 #ifndef __NO_MERSENNE
00075 mt = new unsigned long [MERSENNE_N];
00076 mti = MERSENNE_N + 1;
00077 mersenneMult = 1.0/4294967296.0;
00078 #else
00079 shuffler = new long [NTAB];
00080 #endif
00081 Reset(s);
00082 }
00083
00084 Random::~Random()
00085 {
00086 #ifndef __NO_MERSENNE
00087 delete [] mt;
00088 #else
00089 delete [] shuffler;
00090 #endif
00091 }
00092
00093 void Random::Reset(long s)
00094 {
00095 normSaved = 0;
00096
00097 #ifndef __NO_MERSENNE
00098 InitMersenne(s);
00099 #else
00100
00101 if ((seed = s) < 1)
00102 seed = s == 0 ? 1 : -s;
00103
00104 for (int j=NTAB+7; j>=0; j--)
00105 {
00106 long k = seed / IQ;
00107 seed = IA * (seed - k * IQ) - IR * k;
00108 if (seed < 0) seed += IM;
00109 if (j < NTAB) shuffler[j] = seed;
00110 }
00111 last=shuffler[0];
00112 #endif
00113 }
00114
00115
00116 void Random::InitMersenne(unsigned long s)
00117 {
00118 mt[0]= s & 0xffffffffUL;
00119 for (mti = 1; mti < MERSENNE_N; mti++)
00120 {
00121 mt[mti] = (1812433253UL * (mt[mti-1] ^(mt[mti-1] >> 30)) + mti);
00122
00123
00124
00125
00126
00127 mt[mti] &= 0xffffffffUL;
00128 }
00129 }
00130
00131 int Random::Binary()
00132 {
00133 return Next() > 0.5 ? 1 : 0;
00134 }
00135
00136 #ifndef __NO_MERSENNE
00137
00138 double Random::Next()
00139 {
00140 unsigned long y;
00141
00142
00143 static unsigned long mag01[2]={0x0UL, MATRIX_A};
00144
00145 if (mti >= MERSENNE_N)
00146 {
00147
00148 int kk;
00149
00150
00151 if (mti == MERSENNE_N+1)
00152 InitMersenne(5489UL);
00153
00154 for (kk=0; kk < MERSENNE_N-MERSENNE_M; kk++)
00155 {
00156 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
00157 mt[kk] = mt[kk+MERSENNE_M] ^(y >> 1) ^ mag01[y & 0x1UL];
00158 }
00159
00160 for (; kk < MERSENNE_N-1; kk++)
00161 {
00162 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
00163 mt[kk] = mt[kk+(MERSENNE_M - MERSENNE_N)] ^(y >> 1) ^ mag01[y & 0x1UL];
00164 }
00165
00166 y = (mt[MERSENNE_N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
00167 mt[MERSENNE_N-1] = mt[MERSENNE_M-1] ^(y >> 1) ^ mag01[y & 0x1UL];
00168
00169 mti = 0;
00170 }
00171
00172 y = mt[mti++];
00173
00174
00175 y ^= (y >> 11);
00176 y ^= (y << 7) & 0x9d2c5680UL;
00177 y ^= (y << 15) & 0xefc60000UL;
00178 y ^= (y >> 18);
00179
00180 return (mersenneMult *((double) y + 0.5));
00181 }
00182
00183
00184
00185 unsigned long Random::NextInt()
00186 {
00187 unsigned long y;
00188
00189
00190 static unsigned long mag01[2]={0x0UL, MATRIX_A};
00191
00192 if (mti >= MERSENNE_N)
00193 {
00194
00195 int kk;
00196
00197
00198 if (mti == MERSENNE_N + 1)
00199 InitMersenne(5489UL);
00200
00201 for (kk= 0; kk < MERSENNE_N - MERSENNE_M; kk++)
00202 {
00203 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
00204 mt[kk] = mt[kk+MERSENNE_M] ^(y >> 1) ^ mag01[y & 0x1UL];
00205 }
00206
00207 for (; kk< MERSENNE_N-1; kk++)
00208 {
00209 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
00210 mt[kk] = mt[kk+(MERSENNE_M - MERSENNE_N)] ^(y >> 1) ^ mag01[y & 0x1UL];
00211 }
00212
00213 y = (mt[MERSENNE_N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
00214 mt[MERSENNE_N-1] = mt[MERSENNE_M-1] ^(y >> 1) ^ mag01[y & 0x1UL];
00215
00216 mti = 0;
00217 }
00218
00219 y = mt[mti++];
00220
00221
00222 y ^= (y >> 11);
00223 y ^= (y << 7) & 0x9d2c5680UL;
00224 y ^= (y << 15) & 0xefc60000UL;
00225 y ^= (y >> 18);
00226
00227 return y;
00228 }
00229
00230 #else
00231
00232 double Random::Next()
00233 {
00234
00235
00236 long k = seed / IQ;
00237 seed = IA * (seed - k * IQ) - IR * k;
00238 if (seed < 0) seed += IM;
00239
00240
00241 int j = last/NDIV;
00242
00243
00244 last = shuffler[j];
00245 shuffler[j] = seed;
00246
00247
00248 double temp = AM * last;
00249 if (temp > RNMX) return RNMX;
00250 return temp;
00251 }
00252
00253 unsigned long Random::NextInt()
00254 {
00255
00256
00257 long k = seed / IQ;
00258 seed = IA * (seed - k * IQ) - IR * k;
00259 if (seed < 0) seed += IM;
00260
00261
00262 int j = last/NDIV;
00263
00264
00265 last = shuffler[j];
00266 shuffler[j] = seed;
00267
00268 return last;
00269 }
00270
00271 #endif
00272
00273 double Random::Normal()
00274 {
00275 double v1, v2, fac, rsq;
00276
00277 if (!normSaved)
00278 {
00279 do
00280 {
00281 v1 = 2.0 * Next() - 1.0;
00282 v2 = 2.0 * Next() - 1.0;
00283 rsq = v1*v1 + v2*v2;
00284 }
00285 while (rsq >= 1.0 || rsq == 0.0);
00286
00287 fac = sqrt(-2.0 * log(rsq)/rsq);
00288 normStore = v1 * fac;
00289 normSaved = 1;
00290 return v2 * fac;
00291 }
00292 else
00293 {
00294 normSaved = 0;
00295 return normStore;
00296 }
00297 }
00298
00299 void Random::Choose(int * array, int n, int k)
00300 {
00301 int choices = 1, others = 0;
00302
00303 if (k > n / 2)
00304 {
00305 choices = 0;
00306 others = 1;
00307 k = n - k;
00308 }
00309
00310 for (int i = 0; i < n; i++)
00311 array[i] = others;
00312
00313 while (k > 0)
00314 {
00315 int i = NextInt() % n;
00316
00317 if (array[i] == choices) continue;
00318
00319 array[i] = choices;
00320 k--;
00321 }
00322 }
00323
00324 void Random::Choose(int * array, float * weights, int n, int k)
00325 {
00326 int choices = 1, others = 0;
00327
00328 if (k > n / 2)
00329 {
00330 choices = 0;
00331 others = 1;
00332 k = n - k;
00333 }
00334
00335
00336 float * cumulative = new float [n + 1];
00337
00338 cumulative[0] = 0;
00339 for (int i = 1; i <= n; i++)
00340 cumulative[i] = cumulative[i - 1] + weights[i - 1];
00341
00342 float & sum = cumulative[n], reject = 0.0;
00343
00344 for (int i = 0; i < n; i++)
00345 array[i] = others;
00346
00347 while (k > 0)
00348 {
00349 float weight = Next() * sum;
00350
00351 int hi = n, lo = 0, i = 0;
00352
00353 while (hi >= lo)
00354 {
00355 i = (hi + lo) / 2;
00356
00357 if (cumulative[i + 1] <= weight)
00358 lo = i + 1;
00359 else if (cumulative[i] >= weight)
00360 hi = i - 1;
00361 else break;
00362 }
00363
00364 if (array[i] == choices) continue;
00365
00366 array[i] = choices;
00367 reject += weights[i];
00368
00369
00370
00371 if (reject > sum * 0.50)
00372 {
00373 cumulative[0] = 0;
00374 for (int i = 1; i <= n; i++)
00375 if (array[i] != choices)
00376 cumulative[i] = cumulative[i - 1] + weights[i - 1];
00377 else
00378 cumulative[i] = cumulative[i - 1];
00379
00380 reject = 0.0;
00381 }
00382
00383 k--;
00384 }
00385
00386 delete [] cumulative;
00387 }
00388
00389 Random globalRandom;
00390
00391