Random.cpp

   1
   2 //////////////////////////////////////////////////////////////////////////////
   3 // This file includes code derived from the original Mersenne Twister Code
   4 // by Makoto Matsumoto and Takuji Nishimura
   5 // and is subject to their original copyright notice copied below:
   6 //////////////////////////////////////////////////////////////////////////////
   7
   8 //////////////////////////////////////////////////////////////////////////////
   9 // COPYRIGHT NOTICE FOR MERSENNE TWISTER CODE
  10 // Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
  11 // All rights reserved.
  12 //
  13 // Redistribution and use in source and binary forms, with or without
  14 // modification, are permitted provided that the following conditions
  15 // are met:
  16 //
  17 // 1. Redistributions of source code must retain the above copyright
  18 // notice, this list of conditions and the following disclaimer.
  19 //
  20 // 2. Redistributions in binary form must reproduce the above copyright
  21 // notice, this list of conditions and the following disclaimer in the
  22 // documentation and/or other materials provided with the distribution.
  23 //
  24 // 3. The names of its contributors may not be used to endorse or promote
  25 // products derived from this software without specific prior written
  26 // permission.
  27 //
  28 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  29 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  30 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  31 // A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
  32 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
  33 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
  34 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
  35 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
  36 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
  37 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  38 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  39 //
  40 ///////////////////////////////////////////////////////////////////////////////
  41
  42 #include "Random.h"
  43 #include "MathConstant.h"
  44 #include "Error.h"
  45
  46 #include <math.h>
  47
  48 //Constants used internally by Mersenne random number generator
  49 #define MERSENNE_N 624
  50 #define MERSENNE_M 397
  51
  52 // constant vector a
  53 #define MATRIX_A 0x9908b0dfUL
  54
  55 // most significant w-r bits
  56 #define UPPER_MASK 0x80000000UL
  57
  58 // least significant r bits
  59 #define LOWER_MASK 0x7fffffffUL
  60
  61
  62 // Constants used internally by Park-Miller random generator
  63 #define IA 16807
  64 #define IM 2147483647
  65 #define AM (1.0 / IM)
  66 #define IQ 127773
  67 #define IR 2836
  68 #define NTAB 32
  69 #define NDIV (1+(IM-1)/NTAB)
  70 #define RNMX (1.0-EPS)
  71
  72 Random::Random(long s)
  73    {
  74 #ifndef __NO_MERSENNE
  75    mt = new unsigned long [MERSENNE_N];
  76    mti = MERSENNE_N + 1;
  77    mersenneMult = 1.0/4294967296.0;
  78 #else
  79    shuffler = new long [NTAB];
  80 #endif
  81    Reset(s);
  82    }
  83
  84 Random::~Random()
  85    {
  86 #ifndef __NO_MERSENNE
  87    delete [] mt;
  88 #else
  89    delete [] shuffler;
  90 #endif
  91    }
  92
  93 void Random::Reset(long s)
  94    {
  95    normSaved = 0;
  96
  97 #ifndef __NO_MERSENNE
  98    InitMersenne(s);
  99 #else
 100    // 'Continuous' Random Generator
 101    if ((seed = s) < 1)
 102       seed = s == 0 ? 1 : -s;  // seed == 0 would be disastrous
 103
 104    for (int j=NTAB+7; j>=0; j--)  // Warm up and set shuffle table
 105       {
 106       long k = seed / IQ;
 107       seed = IA * (seed - k * IQ) - IR * k;
 108       if (seed < 0) seed += IM;
 109       if (j < NTAB) shuffler[j] = seed;
 110       }
 111    last=shuffler[0];
 112 #endif
 113    }
 114
 115 // initializes mt[MERSENNE_N] with a seed
 116 void Random::InitMersenne(unsigned long s)
 117    {
 118    mt[0]= s & 0xffffffffUL;
 119    for (mti = 1; mti < MERSENNE_N; mti++)
 120       {
 121       mt[mti] = (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
 122       /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
 123       /* In the previous versions, MSBs of the seed affect */
 124       /* only MSBs of the array mt[].  */
 125       /* 2002/01/09 modified by Makoto Matsumoto */
 126
 127       mt[mti] &= 0xffffffffUL;
 128       }
 129    }
 130
 131 int Random::Binary()
 132    {
 133    return Next() > 0.5 ? 1 : 0;
 134    }
 135
 136 #ifndef __NO_MERSENNE
 137
 138 double Random::Next()
 139    {
 140    unsigned long y;
 141
 142    // mag01[x] = x * MATRIX_A for x=0,1
 143    static unsigned long mag01[2]={0x0UL, MATRIX_A};
 144
 145    if (mti >= MERSENNE_N)
 146       {
 147       /* generate MERSENNE_N words at one time */
 148       int kk;
 149
 150       // If InitMersenne() has not been called, a default initial seed is used
 151       if (mti == MERSENNE_N+1)
 152          InitMersenne(5489UL);
 153
 154       for (kk=0; kk < MERSENNE_N-MERSENNE_M; kk++)
 155          {
 156          y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
 157          mt[kk] = mt[kk+MERSENNE_M] ^ (y >> 1) ^ mag01[y & 0x1UL];
 158          }
 159
 160       for (; kk < MERSENNE_N-1; kk++)
 161          {
 162          y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
 163          mt[kk] = mt[kk+(MERSENNE_M - MERSENNE_N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
 164          }
 165
 166       y = (mt[MERSENNE_N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
 167       mt[MERSENNE_N-1] = mt[MERSENNE_M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
 168
 169       mti = 0;
 170       }
 171
 172    y = mt[mti++];
 173
 174    // Tempering
 175    y ^= (y >> 11);
 176    y ^= (y << 7) & 0x9d2c5680UL;
 177    y ^= (y << 15) & 0xefc60000UL;
 178    y ^= (y >> 18);
 179
 180    return (mersenneMult * ((double) y + 0.5));
 181    }
 182
 183 // Generates a random number on [0,0xffffffff]-interval
 184
 185 unsigned long Random::NextInt()
 186    {
 187    unsigned long y;
 188
 189    // mag01[x] = x * MATRIX_A for x=0,1
 190    static unsigned long mag01[2]={0x0UL, MATRIX_A};
 191
 192    if (mti >= MERSENNE_N)
 193       {
 194       /* generate MERSENNE_N words at one time */
 195       int kk;
 196
 197       // If InitMersenne() has not been called, a default initial seed is used
 198       if (mti == MERSENNE_N + 1)
 199          InitMersenne(5489UL);
 200
 201       for (kk= 0; kk < MERSENNE_N - MERSENNE_M; kk++)
 202          {
 203          y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
 204          mt[kk] = mt[kk+MERSENNE_M] ^ (y >> 1) ^ mag01[y & 0x1UL];
 205          }
 206
 207       for (; kk< MERSENNE_N-1; kk++)
 208          {
 209          y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
 210          mt[kk] = mt[kk+(MERSENNE_M - MERSENNE_N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
 211          }
 212
 213       y = (mt[MERSENNE_N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
 214       mt[MERSENNE_N-1] = mt[MERSENNE_M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
 215
 216       mti = 0;
 217       }
 218
 219    y = mt[mti++];
 220
 221    // Tempering
 222    y ^= (y >> 11);
 223    y ^= (y << 7) & 0x9d2c5680UL;
 224    y ^= (y << 15) & 0xefc60000UL;
 225    y ^= (y >> 18);
 226
 227    return y;
 228    }
 229
 230 #else
 231
 232 double Random::Next()
 233    {
 234    // Compute seed = (IA * seed) % IM without overflows
 235    // by Schrage's method
 236    long k = seed / IQ;
 237    seed = IA * (seed - k * IQ) - IR * k;
 238    if (seed < 0) seed += IM;
 239
 240    // Map to 0..NTAB-1
 241    int j = last/NDIV;
 242
 243    // Output value is shuffler[j], which is in turn replaced by seed
 244    last = shuffler[j];
 245    shuffler[j] = seed;
 246
 247    // Map to 0.0 .. 1.0 excluding endpoints
 248    double temp = AM * last;
 249    if (temp > RNMX) return RNMX;
 250    return temp;
 251    }
 252
 253 unsigned long Random::NextInt()
 254    {
 255    // Compute seed = (IA * seed) % IM without overflows
 256    // by Schrage's method
 257    long k = seed / IQ;
 258    seed = IA * (seed - k * IQ) - IR * k;
 259    if (seed < 0) seed += IM;
 260
 261    // Map to 0..NTAB-1
 262    int j = last/NDIV;
 263
 264    // Output value is shuffler[j], which is in turn replaced by seed
 265    last = shuffler[j];
 266    shuffler[j] = seed;
 267
 268    return last;
 269    }
 270
 271 #endif
 272
 273 double Random::Normal()
 274    {
 275    double v1, v2, fac, rsq;
 276
 277    if (!normSaved)  // Do we need new numbers?
 278       {
 279       do {
 280       v1 = 2.0 * Next() - 1.0;  // Pick two coordinates from
 281       v2 = 2.0 * Next() - 1.0;  // -1 to +1 and check if they
 282       rsq = v1*v1 + v2*v2;  // are in unit circle...
 283       } while (rsq >= 1.0 || rsq == 0.0);
 284
 285       fac = sqrt(-2.0 * log(rsq)/rsq);  // Apply the Box-Muller
 286       normStore = v1 * fac;  // transformation and save
 287       normSaved = 1;  // one deviate for next time
 288       return v2 * fac;
 289       }
 290    else
 291       {
 292       normSaved = 0;
 293       return normStore;
 294       }
 295    }
 296
 297 void Random::Choose(int * array, int n, int k)
 298    {
 299    int choices = 1, others = 0;
 300
 301    if (k > n / 2)
 302       {
 303       choices = 0;
 304       others = 1;
 305       k = n - k;
 306       }
 307
 308    for (int i = 0; i < n; i++)
 309       array[i] = others;
 310
 311    while (k > 0)
 312       {
 313       int i = NextInt() % n;
 314
 315       if (array[i] == choices) continue;
 316
 317       array[i] = choices;
 318       k--;
 319       }
 320    }
 321
 322 void Random::Choose(int * array, float * weights, int n, int k)
 323    {
 324    int choices = 1, others = 0;
 325
 326    if (k > n / 2)
 327       {
 328       choices = 0;
 329       others = 1;
 330       k = n - k;
 331       }
 332
 333    // First calculate cumulative sums of weights ...
 334    float * cumulative = new float [n + 1];
 335
 336    cumulative[0] = 0;
 337    for (int i = 1; i <= n; i++)
 338       cumulative[i] = cumulative[i - 1] + weights[i - 1];
 339
 340    float & sum = cumulative[n], reject = 0.0;
 341
 342    for (int i = 0; i < n; i++)
 343       array[i] = others;
 344
 345    while (k > 0)
 346       {
 347       float weight = (float)(Next() * sum);
 348
 349       int hi = n, lo = 0, i = 0;
 350
 351       while (hi >= lo)
 352          {
 353          i = (hi + lo) / 2;
 354
 355          if (cumulative[i + 1] <= weight)
 356             lo = i + 1;
 357          else if (cumulative[i] >= weight)
 358             hi = i - 1;
 359          else break;
 360          }
 361
 362       if (array[i] == choices) continue;
 363
 364       array[i] = choices;
 365       reject += weights[i];
 366
 367       // After selecting a substantial number of elements, update the cumulative
 368       // distribution -- to ensure that at least half of our samples produce a hit
 369       if (reject > sum * 0.50)
 370          {
 371          cumulative[0] = 0;
 372          for (int i = 1; i <= n; i++)
 373             if (array[i] != choices)
 374                cumulative[i] = cumulative[i - 1] + weights[i - 1];
 375             else
 376                cumulative[i] = cumulative[i - 1];
 377
 378          reject = 0.0;
 379          }
 380
 381       k--;
 382       }
 383    }
 384
 385 Random globalRandom;
 386
 387