libs/basekit/source/RandomGen.c

   1
   2 /*
   3  A C-program for MT19937, with improved initialization 2002/1/26.
   4
   5  This is an optimized version that amortizes the shift/reload cost,
   6  by Eric Landry 2004-03-15.
   7
   8  Before using, initialize the state by using RandomGen_setSeed(seed) or
   9  init_by_array(init_key, key_length).
  10
  11  Copyright (C) 1997--2004, Makoto Matsumoto, Takuji Nishimura, and
  12  Eric Landry; All rights reserved.
  13
  14  Redistribution and use in source and binary forms, with or without
  15  modification, are permitted provided that the following conditions
  16  are met:
  17
  18  1. Redistributions of source code must retain the above copyright
  19  notice, this list of conditions and the following disclaimer.
  20
  21  2. Redistributions in binary form must reproduce the above copyright
  22  notice, this list of conditions and the following disclaimer
  23  in the documentation and/or other materials provided with the
  24  distribution.
  25
  26  3. The names of its contributors may not be used to endorse or
  27  promote products derived from this software without specific
  28  prior written permission.
  29
  30  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  31  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  32  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  33  A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT
  34  OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
  35  SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
  36                                                                           LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  37                                                                           DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
  38  THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  39  (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  40  OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  41
  42  Any feedback is very welcome.
  43  http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
  44  email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
  45
  46  Reference: M. Matsumoto and T. Nishimura, "Mersenne Twister:
  47  A 623-Dimensionally Equidistributed Uniform Pseudo-RandomGen Number
  48  Generator", ACM Transactions on Modeling and Computer Simulation,
  49  Vol. 8, No. 1, January 1998, pp 3--30.
  50  */
  51
  52 #include "Base.h"
  53 #include "RandomGen.h"
  54
  55 /* Period parameters */
  56 #define N RANDOMGEN_N
  57 #define M 397
  58 #define MATRIX_A 0x9908b0dfUL   /* constant vector a */
  59 #define UPPER_MASK 0x80000000UL /* most significant w-r bits */
  60 #define LOWER_MASK 0x7fffffffUL /* least significant r bits */
  61
  62 /*static unsigned long self->mt[N];*/ /* the array for the state vector  */
  63 /*static int self->mti=N+1;*/ /* self->mti==N+1 means self->mt[N] is not initialized */
  64
  65 /* initializes self->mt[N] with a seed */
  66
  67 static void init_genrand(RandomGen *self, unsigned long s)
  68 {
  69         self->mt[0]= s & 0xffffffffUL;
  70         for (self->mti=1; self->mti<N; self->mti ++)
  71         {
  72                 self->mt[self->mti] =
  73                 (1812433253UL * (self->mt[self->mti-1] ^ (self->mt[self->mti-1] >> 30)) + self->mti);
  74                 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
  75                 /* In the previous versions, MSBs of the seed affect   */
  76                 /* only MSBs of the array self->mt[].                        */
  77                 /* 2002/01/09 modified by Makoto Matsumoto             */
  78                 self->mt[self->mti] &= 0xffffffffUL;
  79                 /* for >32 bit machines */
  80         }
  81 }
  82
  83 void RandomGen_setSeed(RandomGen *self, unsigned long seed)
  84 {
  85         init_genrand(self, seed);
  86 }
  87
  88 #include <time.h>
  89
  90 void RandomGen_chooseRandomSeed(RandomGen *self)
  91 {
  92         unsigned long seed = 0;
  93
  94         seed ^= (unsigned long)clock(); // processor time since program start
  95         seed ^= (unsigned long)time(NULL); // seconds since 1970
  96
  97         RandomGen_setSeed(self, seed);
  98 }
  99
 100 /* initialize by an array with array-length */
 101 /* init_key is the array for initializing keys */
 102 /* key_length is its length */
 103 /* slight change for C++, 2004/2/26 */
 104 #ifdef EXTRAS
 105 static void init_by_array(RandomGen *self, unsigned long init_key[],int key_length)
 106 {
 107         int i, j, k;
 108
 109         init_genrand(self, 19650218UL);
 110         i=1; j=0;
 111         k = (N>key_length ? N : key_length);
 112
 113         for (; k; k--) {
 114                 self->mt[i] = (self->mt[i] ^ ((self->mt[i-1] ^ (self->mt[i-1] >> 30)) * 1664525UL))
 115                 + init_key[j] + j; /* non linear */
 116                 self->mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
 117                 i ++; j++;
 118                 if (i>=N) { self->mt[0] = self->mt[N-1]; i=1; }
 119                 if (j>=key_length) j=0;
 120         }
 121
 122         for (k=N-1; k; k--) {
 123                 self->mt[i] = (self->mt[i] ^ ((self->mt[i-1] ^ (self->mt[i-1] >> 30)) * 1566083941UL))
 124                 - i; /* non linear */
 125                 self->mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
 126                 i ++;
 127                 if (i>=N) { self->mt[0] = self->mt[N-1]; i=1; }
 128         }
 129
 130         self->mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
 131 }
 132 #endif
 133
 134 /* generates a random number on [0,0xffffffff]-interval */
 135 static unsigned long genrand_int32(RandomGen *self)
 136 {
 137         unsigned long y;
 138         static unsigned long mag01[2]={0x0UL, MATRIX_A};
 139         /* mag01[x] = x * MATRIX_A  for x=0,1 */
 140
 141         if (self->mti >= N) { /* generate N words at one time */
 142                 int kk;
 143
 144                 if (self->mti == N+1)   /* if init_genrand() has not been called, */
 145                         init_genrand(self, 5489UL); /* a default initial seed is used */
 146
 147                 for (kk=0;kk<N-M;kk++) {
 148                         y = (self->mt[kk]&UPPER_MASK)|(self->mt[kk+1]&LOWER_MASK);
 149                         self->mt[kk] = self->mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];
 150                 }
 151                 for (;kk<N-1;kk++) {
 152                         y = (self->mt[kk]&UPPER_MASK)|(self->mt[kk+1]&LOWER_MASK);
 153                         self->mt[kk] = self->mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
 154                 }
 155                 y = (self->mt[N-1]&UPPER_MASK)|(self->mt[0]&LOWER_MASK);
 156                 self->mt[N-1] = self->mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
 157
 158                 self->mti = 0;
 159         }
 160
 161         y = self->mt[self->mti ++];
 162
 163         // Tempering
 164         y ^= (y >> 11);
 165         y ^= (y << 7) & 0x9d2c5680UL;
 166         y ^= (y << 15) & 0xefc60000UL;
 167         y ^= (y >> 18);
 168
 169         return y;
 170 }
 171
 172 /* generates a random number on [0,0x7fffffff]-interval */
 173 #ifdef EXTRAS
 174 static long genrand_int31(RandomGen *self)
 175 {
 176         return (long)(genrand_int32(self)>>1);
 177 }
 178
 179 /* generates a random number on [0,1]-real-interval */
 180 static double genrand_real1(RandomGen *self)
 181 {
 182         return genrand_int32(self)*(1.0/4294967295.0);
 183         /* divided by 2^32-1 */
 184 }
 185 #endif
 186
 187 /* generates a random number on [0,1)-real-interval */
 188 static double genrand_real2(RandomGen *self)
 189 {
 190         return genrand_int32(self)*(1.0/4294967296.0);
 191         /* divided by 2^32 */
 192 }
 193
 194 /* generates a random number on (0,1)-real-interval */
 195 #ifdef EXTRAS
 196 static double genrand_real3(RandomGen *self)
 197 {
 198         return (((double)genrand_int32(self)) + 0.5)*(1.0/4294967296.0);
 199         /* divided by 2^32 */
 200 }
 201
 202 /* generates a random number on [0,1) with 53-bit resolution*/
 203 static double genrand_res53(RandomGen *self)
 204 {
 205         unsigned long a=genrand_int32(self)>>5, b=genrand_int32(self)>>6;
 206         return(a*67108864.0+b)*(1.0/9007199254740992.0);
 207 }
 208 #endif
 209
 210 /* These real versions are due to Isaku Wada, 2002/01/09 added */
 211
 212 /*
 213  int main(void)
 214  {
 215          int i;
 216          unsigned long init[4]={0x123, 0x234, 0x345, 0x456}, length=4;
 217          init_by_array(init, length);
 218          printf("1000 outputs of genrand_int32(self)\n");
 219          for (i=0; i<1000; i ++) {
 220                  printf("%10lu ", genrand_int32(self));
 221                  if (i%5==4) printf("\n");
 222          }
 223          printf("\n1000 outputs of genrand_real2()\n");
 224          for (i=0; i<1000; i ++) {
 225                  printf("%10.8f ", genrand_real2());
 226                  if (i%5==4) printf("\n");
 227          }
 228          return 0;
 229  }
 230  */
 231
 232 /* --------------------------------------------------------
 233    This stuff added by Steve Dekorte, 2004 07 04
 234 */
 235
 236 #include <time.h>
 237 #include <stdlib.h>
 238
 239 RandomGen *RandomGen_new(void)
 240 {
 241         RandomGen *self = (RandomGen *)io_calloc(1, sizeof(RandomGen));
 242         unsigned long t1 = (unsigned long)time(NULL);
 243         unsigned long t2 = clock();
 244         self->mti = RANDOMGEN_N + 1;
 245         init_genrand(self, t1 + t2);
 246         self->y2 = 0;
 247         return self;
 248 }
 249
 250 void RandomGen_free(RandomGen *self)
 251 {
 252         io_free(self);
 253 }
 254
 255 double RandomGen_randomDouble(RandomGen *self)
 256 {
 257         return genrand_real2(self);
 258 }
 259
 260 int RandomGen_randomInt(RandomGen *self)
 261 {
 262         return genrand_int32(self);
 263 }
 264
 265 #include <math.h>
 266 #ifndef M_PI_2 // some windows'es don't define this
 267 # define M_PI_2 1.57079632679489661923
 268 #endif
 269
 270 double RandomGen_gaussian(RandomGen *self, double m, double s)
 271 {
 272         // http://www.taygeta.com/random/gaussian.html
 273         double x1, x2, w, y1, y2;
 274
 275         do {
 276                 x1 = 2.0 * RandomGen_randomDouble(self) - 1.0;
 277                 x2 = 2.0 * RandomGen_randomDouble(self) - 1.0;
 278                 w = x1 * x1 + x2 * x2;
 279         } while ( w >= 1.0 );
 280
 281         w = sqrt( (-2.0 * log( w ) ) / w );
 282         y1 = x1 * w;
 283         y2 = x2 * w;
 284
 285         // The following code resulted in a lot of nans being returned. The
 286         // following code *should* also be slower.
 287         /*
 288         double x1 = 2.0 * RandomGen_randomDouble(self) - 1.0;
 289         double x2 = 2.0 * RandomGen_randomDouble(self) - 1.0;
 290         double y1 = sqrt( - 2.0 * log(x1) ) * cos( M_PI_2 * x2 );
 291         */
 292
 293         return ( m + y1 * s );
 294 }
 295
 296 /*
 297 double RandomGen_gaussian(RandomGen *self, double m, double s)
 298 {
 299         // mean m, standard deviation s
 300         double x1, x2, w, y1;
 301
 302         if (self->use_last) // use value from previous call
 303         {
 304                 y1 = self->y2;
 305                 self->use_last = 0;
 306         }
 307         else
 308         {
 309                 do {
 310                         x1 = 2.0 * RandomGen_randomDouble(self) - 1.0;
 311                         x2 = 2.0 * RandomGen_randomDouble(self) - 1.0;
 312                         w = x1 * x1 + x2 * x2;
 313                 } while ( w >= 1.0 );
 314
 315                 w = sqrt( (-2.0 * log( w ) ) / w );
 316                 y1 = x1 * w;
 317                 self->y2 = x2 * w;
 318                 self->use_last = 1;
 319         }
 320
 321         return ( m + y1 * s );
 322 }
 323 */
 324