random.c

   1 /**********************************************************************
   2
   3   random.c -
   4
   5   $Author$
   6   created at: Fri Dec 24 16:39:21 JST 1993
   7
   8   Copyright (C) 1993-2007 Yukihiro Matsumoto
   9
  10 **********************************************************************/
  11
  12 /*
  13 This is based on trimmed version of MT19937.  To get the original version,
  14 contact <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html>.
  15
  16 The original copyright notice follows.
  17
  18    A C-program for MT19937, with initialization improved 2002/2/10.
  19    Coded by Takuji Nishimura and Makoto Matsumoto.
  20    This is a faster version by taking Shawn Cokus's optimization,
  21    Matthe Bellew's simplification, Isaku Wada's real version.
  22
  23    Before using, initialize the state by using init_genrand(seed)
  24    or init_by_array(init_key, key_length).
  25
  26    Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
  27    All rights reserved.
  28
  29    Redistribution and use in source and binary forms, with or without
  30    modification, are permitted provided that the following conditions
  31    are met:
  32
  33      1. Redistributions of source code must retain the above copyright
  34         notice, this list of conditions and the following disclaimer.
  35
  36      2. Redistributions in binary form must reproduce the above copyright
  37         notice, this list of conditions and the following disclaimer in the
  38         documentation and/or other materials provided with the distribution.
  39
  40      3. The names of its contributors may not be used to endorse or promote
  41         products derived from this software without specific prior written
  42         permission.
  43
  44    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  45    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  46    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  47    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
  48    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
  49    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
  50    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
  51    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
  52    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
  53    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  54    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  55
  56
  57    Any feedback is very welcome.
  58    http://www.math.keio.ac.jp/matumoto/emt.html
  59    email: matumoto@math.keio.ac.jp
  60 */
  61
  62 /* Period parameters */
  63 #define N 624
  64 #define M 397
  65 #define MATRIX_A 0x9908b0dfUL   /* constant vector a */
  66 #define UMASK 0x80000000UL /* most significant w-r bits */
  67 #define LMASK 0x7fffffffUL /* least significant r bits */
  68 #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) )
  69 #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL))
  70
  71 static unsigned long state[N]; /* the array for the state vector  */
  72 static int left = 1;
  73 static int initf = 0;
  74 static unsigned long *next;
  75
  76 /* initializes state[N] with a seed */
  77 static void
  78 init_genrand(unsigned long s)
  79 {
  80     int j;
  81     state[0]= s & 0xffffffffUL;
  82     for (j=1; j<N; j++) {
  83         state[j] = (1812433253UL * (state[j-1] ^ (state[j-1] >> 30)) + j);
  84         /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
  85         /* In the previous versions, MSBs of the seed affect   */
  86         /* only MSBs of the array state[].                        */
  87         /* 2002/01/09 modified by Makoto Matsumoto             */
  88         state[j] &= 0xffffffffUL;  /* for >32 bit machines */
  89     }
  90     left = 1; initf = 1;
  91 }
  92
  93 /* initialize by an array with array-length */
  94 /* init_key is the array for initializing keys */
  95 /* key_length is its length */
  96 /* slight change for C++, 2004/2/26 */
  97 static void
  98 init_by_array(unsigned long init_key[], int key_length)
  99 {
 100     int i, j, k;
 101     init_genrand(19650218UL);
 102     i=1; j=0;
 103     k = (N>key_length ? N : key_length);
 104     for (; k; k--) {
 105         state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1664525UL))
 106           + init_key[j] + j; /* non linear */
 107         state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
 108         i++; j++;
 109         if (i>=N) { state[0] = state[N-1]; i=1; }
 110         if (j>=key_length) j=0;
 111     }
 112     for (k=N-1; k; k--) {
 113         state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL))
 114           - i; /* non linear */
 115         state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
 116         i++;
 117         if (i>=N) { state[0] = state[N-1]; i=1; }
 118     }
 119
 120     state[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
 121     left = 1; initf = 1;
 122 }
 123
 124 static void
 125 next_state(void)
 126 {
 127     unsigned long *p=state;
 128     int j;
 129
 130     /* if init_genrand() has not been called, */
 131     /* a default initial seed is used         */
 132     if (initf==0) init_genrand(5489UL);
 133
 134     left = N;
 135     next = state;
 136
 137     for (j=N-M+1; --j; p++)
 138         *p = p[M] ^ TWIST(p[0], p[1]);
 139
 140     for (j=M; --j; p++)
 141         *p = p[M-N] ^ TWIST(p[0], p[1]);
 142
 143     *p = p[M-N] ^ TWIST(p[0], state[0]);
 144 }
 145
 146 /* generates a random number on [0,0xffffffff]-interval */
 147 static unsigned long
 148 genrand_int32(void)
 149 {
 150     unsigned long y;
 151
 152     if (--left == 0) next_state();
 153     y = *next++;
 154
 155     /* Tempering */
 156     y ^= (y >> 11);
 157     y ^= (y << 7) & 0x9d2c5680UL;
 158     y ^= (y << 15) & 0xefc60000UL;
 159     y ^= (y >> 18);
 160
 161     return y;
 162 }
 163
 164 /* generates a random number on [0,1) with 53-bit resolution*/
 165 static double
 166 genrand_real(void)
 167 {
 168     unsigned long a=genrand_int32()>>5, b=genrand_int32()>>6;
 169     return(a*67108864.0+b)*(1.0/9007199254740992.0);
 170 }
 171 /* These real versions are due to Isaku Wada, 2002/01/09 added */
 172
 173 #undef N
 174 #undef M
 175
 176 /* These real versions are due to Isaku Wada, 2002/01/09 added */
 177
 178 #include "ruby/ruby.h"
 179
 180 #ifdef HAVE_UNISTD_H
 181 #include <unistd.h>
 182 #endif
 183 #include <time.h>
 184 #include <sys/types.h>
 185 #include <sys/stat.h>
 186 #ifdef HAVE_FCNTL_H
 187 #include <fcntl.h>
 188 #endif
 189
 190 unsigned long
 191 rb_genrand_int32(void)
 192 {
 193     return genrand_int32();
 194 }
 195
 196 double
 197 rb_genrand_real(void)
 198 {
 199     return genrand_real();
 200 }
 201
 202 static int seed_initialized = 0;
 203 static VALUE saved_seed = INT2FIX(0);
 204
 205 static VALUE
 206 rand_init(VALUE vseed)
 207 {
 208     volatile VALUE seed;
 209     VALUE old;
 210     long len;
 211     unsigned long *buf;
 212
 213     seed = rb_to_int(vseed);
 214     switch (TYPE(seed)) {
 215       case T_FIXNUM:
 216           len = sizeof(VALUE);
 217           break;
 218       case T_BIGNUM:
 219           len = RBIGNUM_LEN(seed) * SIZEOF_BDIGITS;
 220           if (len == 0)
 221               len = 4;
 222           break;
 223       default:
 224           rb_raise(rb_eTypeError, "failed to convert %s into Integer",
 225                    rb_obj_classname(vseed));
 226     }
 227     len = (len + 3) / 4; /* number of 32bit words */
 228     buf = ALLOC_N(unsigned long, len); /* allocate longs for init_by_array */
 229     memset(buf, 0, len * sizeof(long));
 230     if (FIXNUM_P(seed)) {
 231         buf[0] = FIX2ULONG(seed) & 0xffffffff;
 232 #if SIZEOF_LONG > 4
 233         buf[1] = FIX2ULONG(seed) >> 32;
 234 #endif
 235     }
 236     else {
 237         int i, j;
 238         for (i = RBIGNUM_LEN(seed)-1; 0 <= i; i--) {
 239             j = i * SIZEOF_BDIGITS / 4;
 240 #if SIZEOF_BDIGITS < 4
 241             buf[j] <<= SIZEOF_BDIGITS * 8;
 242 #endif
 243             buf[j] |= RBIGNUM_DIGITS(seed)[i];
 244         }
 245     }
 246     while (1 < len && buf[len-1] == 0) {
 247         len--;
 248     }
 249     if (len <= 1) {
 250         init_genrand(buf[0]);
 251     }
 252     else {
 253         if (buf[len-1] == 1) /* remove leading-zero-guard */
 254             len--;
 255         init_by_array(buf, len);
 256     }
 257     old = saved_seed;
 258     saved_seed = seed;
 259     xfree(buf);
 260     return old;
 261 }
 262
 263 #define DEFAULT_SEED_LEN (4 * sizeof(long))
 264
 265 static void
 266 fill_random_seed(char *ptr)
 267 {
 268     static int n = 0;
 269     unsigned long *seed;
 270     struct timeval tv;
 271     int fd;
 272     struct stat statbuf;
 273     char *buf = (char*)ptr;
 274
 275     seed = (unsigned long *)buf;
 276
 277     memset(buf, 0, DEFAULT_SEED_LEN);
 278
 279 #ifdef S_ISCHR
 280     if ((fd = open("/dev/urandom", O_RDONLY
 281 #ifdef O_NONBLOCK
 282             |O_NONBLOCK
 283 #endif
 284 #ifdef O_NOCTTY
 285             |O_NOCTTY
 286 #endif
 287 #ifdef O_NOFOLLOW
 288             |O_NOFOLLOW
 289 #endif
 290             )) >= 0) {
 291         if (fstat(fd, &statbuf) == 0 && S_ISCHR(statbuf.st_mode)) {
 292             read(fd, seed, DEFAULT_SEED_LEN);
 293         }
 294         close(fd);
 295     }
 296 #endif
 297
 298     gettimeofday(&tv, 0);
 299     seed[0] ^= tv.tv_usec;
 300     seed[1] ^= tv.tv_sec;
 301     seed[2] ^= getpid() ^ (n++ << 16);
 302     seed[3] ^= (unsigned long)&seed;
 303 }
 304
 305 static VALUE
 306 make_seed_value(char *ptr)
 307 {
 308     BDIGIT *digits;
 309     NEWOBJ(big, struct RBignum);
 310     OBJSETUP(big, rb_cBignum, T_BIGNUM);
 311
 312     RBIGNUM_SET_SIGN(big, 1);
 313     rb_big_resize((VALUE)big, DEFAULT_SEED_LEN / SIZEOF_BDIGITS + 1);
 314     digits = RBIGNUM_DIGITS(big);
 315
 316     MEMCPY((char *)RBIGNUM_DIGITS(big), ptr, char, DEFAULT_SEED_LEN);
 317
 318     /* set leading-zero-guard if need. */
 319     digits[RBIGNUM_LEN(big)-1] = digits[RBIGNUM_LEN(big)-2] <= 1 ? 1 : 0;
 320
 321     return rb_big_norm((VALUE)big);
 322 }
 323
 324 static VALUE
 325 random_seed(void)
 326 {
 327     char buf[DEFAULT_SEED_LEN];
 328     fill_random_seed(buf);
 329     return make_seed_value(buf);
 330 }
 331
 332 /*
 333  *  call-seq:
 334  *     srand(number=0)    => old_seed
 335  *
 336  *  Seeds the pseudorandom number generator to the value of
 337  *  <i>number</i>.<code>to_i.abs</code>. If <i>number</i> is omitted
 338  *  or zero, seeds the generator using a combination of the time, the
 339  *  process id, and a sequence number. (This is also the behavior if
 340  *  <code>Kernel::rand</code> is called without previously calling
 341  *  <code>srand</code>, but without the sequence.) By setting the seed
 342  *  to a known value, scripts can be made deterministic during testing.
 343  *  The previous seed value is returned. Also see <code>Kernel::rand</code>.
 344  */
 345
 346 static VALUE
 347 rb_f_srand(int argc, VALUE *argv, VALUE obj)
 348 {
 349     VALUE seed, old;
 350
 351     rb_secure(4);
 352     if (argc == 0) {
 353         seed = random_seed();
 354     }
 355     else {
 356         rb_scan_args(argc, argv, "01", &seed);
 357     }
 358     old = rand_init(seed);
 359
 360     return old;
 361 }
 362
 363 static unsigned long
 364 make_mask(unsigned long x)
 365 {
 366     x = x | x >> 1;
 367     x = x | x >> 2;
 368     x = x | x >> 4;
 369     x = x | x >> 8;
 370     x = x | x >> 16;
 371 #if 4 < SIZEOF_LONG
 372     x = x | x >> 32;
 373 #endif
 374     return x;
 375 }
 376
 377 static unsigned long
 378 limited_rand(unsigned long limit)
 379 {
 380     unsigned long mask = make_mask(limit);
 381     int i;
 382     unsigned long val;
 383
 384   retry:
 385     val = 0;
 386     for (i = SIZEOF_LONG/4-1; 0 <= i; i--) {
 387         if (mask >> (i * 32)) {
 388             val |= genrand_int32() << (i * 32);
 389             val &= mask;
 390             if (limit < val)
 391                 goto retry;
 392         }
 393     }
 394     return val;
 395 }
 396
 397 static VALUE
 398 limited_big_rand(struct RBignum *limit)
 399 {
 400     unsigned long mask, lim, rnd;
 401     struct RBignum *val;
 402     int i, len, boundary;
 403
 404     len = (RBIGNUM_LEN(limit) * SIZEOF_BDIGITS + 3) / 4;
 405     val = (struct RBignum *)rb_big_clone((VALUE)limit);
 406     RBIGNUM_SET_SIGN(val, 1);
 407 #if SIZEOF_BDIGITS == 2
 408 # define BIG_GET32(big,i) \
 409     (RBIGNUM_DIGITS(big)[(i)*2] | \
 410      ((i)*2+1 < RBIGNUM_LEN(big) ? \
 411       (RBIGNUM_DIGITS(big)[(i)*2+1] << 16) : \
 412       0))
 413 # define BIG_SET32(big,i,d) \
 414     ((RBIGNUM_DIGITS(big)[(i)*2] = (d) & 0xffff), \
 415      ((i)*2+1 < RBIGNUM_LEN(big) ? \
 416       (RBIGNUM_DIGITS(big)[(i)*2+1] = (d) >> 16) : \
 417       0))
 418 #else
 419     /* SIZEOF_BDIGITS == 4 */
 420 # define BIG_GET32(big,i) (RBIGNUM_DIGITS(big)[i])
 421 # define BIG_SET32(big,i,d) (RBIGNUM_DIGITS(big)[i] = (d))
 422 #endif
 423   retry:
 424     mask = 0;
 425     boundary = 1;
 426     for (i = len-1; 0 <= i; i--) {
 427         lim = BIG_GET32(limit, i);
 428         mask = mask ? 0xffffffff : make_mask(lim);
 429         if (mask) {
 430             rnd = genrand_int32() & mask;
 431             if (boundary) {
 432                 if (lim < rnd)
 433                     goto retry;
 434                 if (rnd < lim)
 435                     boundary = 0;
 436             }
 437         }
 438         else {
 439             rnd = 0;
 440         }
 441         BIG_SET32(val, i, rnd);
 442     }
 443     return rb_big_norm((VALUE)val);
 444 }
 445
 446 /*
 447  *  call-seq:
 448  *     rand(max=0)    => number
 449  *
 450  *  Converts <i>max</i> to an integer using max1 =
 451  *  max<code>.to_i.abs</code>. If the result is zero, returns a
 452  *  pseudorandom floating point number greater than or equal to 0.0 and
 453  *  less than 1.0. Otherwise, returns a pseudorandom integer greater
 454  *  than or equal to zero and less than max1. <code>Kernel::srand</code>
 455  *  may be used to ensure repeatable sequences of random numbers between
 456  *  different runs of the program. Ruby currently uses a modified
 457  *  Mersenne Twister with a period of 2**19937-1.
 458  *
 459  *     srand 1234                 #=> 0
 460  *     [ rand,  rand ]            #=> [0.191519450163469, 0.49766366626136]
 461  *     [ rand(10), rand(1000) ]   #=> [6, 817]
 462  *     srand 1234                 #=> 1234
 463  *     [ rand,  rand ]            #=> [0.191519450163469, 0.49766366626136]
 464  */
 465
 466 static VALUE
 467 rb_f_rand(int argc, VALUE *argv, VALUE obj)
 468 {
 469     VALUE vmax;
 470     long val, max;
 471
 472     rb_scan_args(argc, argv, "01", &vmax);
 473     if (!seed_initialized) {
 474        rand_init(random_seed());
 475     }
 476     switch (TYPE(vmax)) {
 477       case T_FLOAT:
 478         if (RFLOAT_VALUE(vmax) <= LONG_MAX && RFLOAT_VALUE(vmax) >= LONG_MIN) {
 479             max = (long)RFLOAT_VALUE(vmax);
 480             break;
 481         }
 482         if (RFLOAT_VALUE(vmax) < 0)
 483             vmax = rb_dbl2big(-RFLOAT_VALUE(vmax));
 484         else
 485             vmax = rb_dbl2big(RFLOAT_VALUE(vmax));
 486         /* fall through */
 487       case T_BIGNUM:
 488       bignum:
 489         {
 490             struct RBignum *limit = (struct RBignum *)vmax;
 491             if (!RBIGNUM_SIGN(limit)) {
 492                 limit = (struct RBignum *)rb_big_clone(vmax);
 493                 RBIGNUM_SET_SIGN(limit, 1);
 494             }
 495             limit = (struct RBignum *)rb_big_minus((VALUE)limit, INT2FIX(1));
 496             if (FIXNUM_P((VALUE)limit)) {
 497                 if (FIX2LONG((VALUE)limit) == -1)
 498                     return DOUBLE2NUM(genrand_real());
 499                 return LONG2NUM(limited_rand(FIX2LONG((VALUE)limit)));
 500             }
 501             return limited_big_rand(limit);
 502         }
 503       case T_NIL:
 504         max = 0;
 505         break;
 506       default:
 507         vmax = rb_Integer(vmax);
 508         if (TYPE(vmax) == T_BIGNUM) goto bignum;
 509       case T_FIXNUM:
 510         max = FIX2LONG(vmax);
 511         break;
 512     }
 513
 514     if (max == 0) {
 515         return DOUBLE2NUM(genrand_real());
 516     }
 517     if (max < 0) max = -max;
 518     val = limited_rand(max-1);
 519     return LONG2NUM(val);
 520 }
 521
 522 static char initial_seed[DEFAULT_SEED_LEN];
 523
 524 void
 525 Init_RandomSeed(void)
 526 {
 527     fill_random_seed(initial_seed);
 528     init_by_array((unsigned long*)initial_seed, DEFAULT_SEED_LEN/sizeof(unsigned long));
 529     seed_initialized = 1;
 530 }
 531
 532 static void
 533 Init_RandomSeed2(void)
 534 {
 535     saved_seed = make_seed_value(initial_seed);
 536     memset(initial_seed, 0, DEFAULT_SEED_LEN);
 537 }
 538
 539 void
 540 rb_reset_random_seed(void)
 541 {
 542     seed_initialized = 0;
 543     saved_seed = INT2FIX(0);
 544 }
 545
 546 void
 547 Init_Random(void)
 548 {
 549     Init_RandomSeed2();
 550     rb_define_global_function("srand", rb_f_srand, -1);
 551     rb_define_global_function("rand", rb_f_rand, -1);
 552     rb_global_variable(&saved_seed);
 553 }