60#include "MathConstant.h"
70#define MATRIX_A 0x9908b0dfUL
73#define UPPER_MASK 0x80000000UL
76#define LOWER_MASK 0x7fffffffUL
86#define NDIV (1+(IM-1)/NTAB)
92 mt =
new unsigned long [MERSENNE_N];
94 mersenneMult = 1.0/4294967296.0;
96 shuffler =
new long [NTAB];
110void Random::Reset(
long s)
119 seed = s == 0 ? 1 : -s;
121 for (
int j=NTAB+7; j>=0; j--)
124 seed = IA * (seed - k * IQ) - IR * k;
125 if (seed < 0) seed += IM;
126 if (j < NTAB) shuffler[j] = seed;
133void Random::InitMersenne(
unsigned long s)
135 mt[0]= s & 0xffffffffUL;
136 for (mti = 1; mti < MERSENNE_N; mti++)
138 mt[mti] = (1812433253UL * (mt[mti-1] ^(mt[mti-1] >> 30)) + mti);
144 mt[mti] &= 0xffffffffUL;
150 return Next() > 0.5 ? 1 : 0;
160 static unsigned long mag01[2]={0x0UL, MATRIX_A};
162 if (mti >= MERSENNE_N)
168 if (mti == MERSENNE_N+1)
169 InitMersenne(5489UL);
171 for (kk=0; kk < MERSENNE_N-MERSENNE_M; kk++)
173 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
174 mt[kk] = mt[kk+MERSENNE_M] ^(y >> 1) ^ mag01[y & 0x1UL];
177 for (; kk < MERSENNE_N-1; kk++)
179 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
180 mt[kk] = mt[kk+(MERSENNE_M - MERSENNE_N)] ^(y >> 1) ^ mag01[y & 0x1UL];
183 y = (mt[MERSENNE_N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
184 mt[MERSENNE_N-1] = mt[MERSENNE_M-1] ^(y >> 1) ^ mag01[y & 0x1UL];
193 y ^= (y << 7) & 0x9d2c5680UL;
194 y ^= (y << 15) & 0xefc60000UL;
197 return (mersenneMult *((
double) y + 0.5));
202unsigned long Random::NextInt()
207 static unsigned long mag01[2]={0x0UL, MATRIX_A};
209 if (mti >= MERSENNE_N)
215 if (mti == MERSENNE_N + 1)
216 InitMersenne(5489UL);
218 for (kk= 0; kk < MERSENNE_N - MERSENNE_M; kk++)
220 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
221 mt[kk] = mt[kk+MERSENNE_M] ^(y >> 1) ^ mag01[y & 0x1UL];
224 for (; kk< MERSENNE_N-1; kk++)
226 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
227 mt[kk] = mt[kk+(MERSENNE_M - MERSENNE_N)] ^(y >> 1) ^ mag01[y & 0x1UL];
230 y = (mt[MERSENNE_N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
231 mt[MERSENNE_N-1] = mt[MERSENNE_M-1] ^(y >> 1) ^ mag01[y & 0x1UL];
240 y ^= (y << 7) & 0x9d2c5680UL;
241 y ^= (y << 15) & 0xefc60000UL;
254 seed = IA * (seed - k * IQ) - IR * k;
255 if (seed < 0) seed += IM;
265 double temp = AM * last;
266 if (temp > RNMX)
return RNMX;
270unsigned long Random::NextInt()
275 seed = IA * (seed - k * IQ) - IR * k;
276 if (seed < 0) seed += IM;
290double Random::Normal()
292 double v1, v2, fac, rsq;
298 v1 = 2.0 * Next() - 1.0;
299 v2 = 2.0 * Next() - 1.0;
302 while (rsq >= 1.0 || rsq == 0.0);
304 fac = sqrt(-2.0 * log(rsq)/rsq);
305 normStore = v1 * fac;
316void Random::Choose(
int * array,
int n,
int k)
318 int choices = 1, others = 0;
327 for (
int i = 0; i < n; i++)
332 int i = NextInt() % n;
334 if (array[i] == choices)
continue;
341void Random::Choose(
int * array,
float * weights,
int n,
int k)
343 int choices = 1, others = 0;
353 float * cumulative =
new float [n + 1];
356 for (
int i = 1; i <= n; i++)
357 cumulative[i] = cumulative[i - 1] + weights[i - 1];
359 float & sum = cumulative[n], reject = 0.0;
361 for (
int i = 0; i < n; i++)
366 float weight = Next() * sum;
368 int hi = n, lo = 0, i = 0;
374 if (cumulative[i + 1] <= weight)
376 else if (cumulative[i] >= weight)
381 if (array[i] == choices)
continue;
384 reject += weights[i];
388 if (reject > sum * 0.50)
391 for (
int i = 1; i <= n; i++)
392 if (array[i] != choices)
393 cumulative[i] = cumulative[i - 1] + weights[i - 1];
395 cumulative[i] = cumulative[i - 1];
404 delete [] cumulative;