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 VALUE saved_seed = INT2FIX(0);
 203
 204 static VALUE
 205 rand_init(VALUE vseed)
 206 {
 207     volatile VALUE seed;
 208     VALUE old;
 209     long len;
 210     unsigned long *buf;
 211
 212     seed = rb_to_int(vseed);
 213     switch (TYPE(seed)) {
 214       case T_FIXNUM:
 215           len = sizeof(VALUE);
 216           break;
 217       case T_BIGNUM:
 218           len = RBIGNUM_LEN(seed) * SIZEOF_BDIGITS;
 219           if (len == 0)
 220               len = 4;
 221           break;
 222       default:
 223           rb_raise(rb_eTypeError, "failed to convert %s into Integer",
 224                    rb_obj_classname(vseed));
 225     }
 226     len = (len + 3) / 4; /* number of 32bit words */
 227     buf = ALLOC_N(unsigned long, len); /* allocate longs for init_by_array */
 228     memset(buf, 0, len * sizeof(long));
 229     if (FIXNUM_P(seed)) {
 230         buf[0] = FIX2ULONG(seed) & 0xffffffff;
 231 #if SIZEOF_LONG > 4
 232         buf[1] = FIX2ULONG(seed) >> 32;
 233 #endif
 234     }
 235     else {
 236         int i, j;
 237         for (i = RBIGNUM_LEN(seed)-1; 0 <= i; i--) {
 238             j = i * SIZEOF_BDIGITS / 4;
 239 #if SIZEOF_BDIGITS < 4
 240             buf[j] <<= SIZEOF_BDIGITS * 8;
 241 #endif
 242             buf[j] |= RBIGNUM_DIGITS(seed)[i];
 243         }
 244     }
 245     while (1 < len && buf[len-1] == 0) {
 246         len--;
 247     }
 248     if (len <= 1) {
 249         init_genrand(buf[0]);
 250     }
 251     else {
 252         if (buf[len-1] == 1) /* remove leading-zero-guard */
 253             len--;
 254         init_by_array(buf, len);
 255     }
 256     old = saved_seed;
 257     saved_seed = seed;
 258     free(buf);
 259     return old;
 260 }
 261
 262 static VALUE
 263 random_seed(void)
 264 {
 265     static int n = 0;
 266     struct timeval tv;
 267     int fd;
 268     struct stat statbuf;
 269
 270     int seed_len;
 271     BDIGIT *digits;
 272     unsigned long *seed;
 273     NEWOBJ(big, struct RBignum);
 274     OBJSETUP(big, rb_cBignum, T_BIGNUM);
 275
 276     seed_len = 4 * sizeof(long);
 277     RBIGNUM_SET_SIGN(big, 1);
 278     rb_big_resize((VALUE)big, seed_len / SIZEOF_BDIGITS + 1);
 279     digits = RBIGNUM_DIGITS(big);
 280     seed = (unsigned long *)RBIGNUM_DIGITS(big);
 281
 282     memset(digits, 0, RBIGNUM_LEN(big) * SIZEOF_BDIGITS);
 283
 284 #ifdef S_ISCHR
 285     if ((fd = open("/dev/urandom", O_RDONLY
 286 #ifdef O_NONBLOCK
 287             |O_NONBLOCK
 288 #endif
 289 #ifdef O_NOCTTY
 290             |O_NOCTTY
 291 #endif
 292 #ifdef O_NOFOLLOW
 293             |O_NOFOLLOW
 294 #endif
 295             )) >= 0) {
 296         if (fstat(fd, &statbuf) == 0 && S_ISCHR(statbuf.st_mode)) {
 297             read(fd, seed, seed_len);
 298         }
 299         close(fd);
 300     }
 301 #endif
 302
 303     gettimeofday(&tv, 0);
 304     seed[0] ^= tv.tv_usec;
 305     seed[1] ^= tv.tv_sec;
 306     seed[2] ^= getpid() ^ (n++ << 16);
 307     seed[3] ^= (unsigned long)&seed;
 308
 309     /* set leading-zero-guard if need. */
 310     digits[RBIGNUM_LEN(big)-1] = digits[RBIGNUM_LEN(big)-2] <= 1 ? 1 : 0;
 311
 312     return rb_big_norm((VALUE)big);
 313 }
 314
 315 /*
 316  *  call-seq:
 317  *     srand(number=0)    => old_seed
 318  *
 319  *  Seeds the pseudorandom number generator to the value of
 320  *  <i>number</i>.<code>to_i.abs</code>. If <i>number</i> is omitted
 321  *  or zero, seeds the generator using a combination of the time, the
 322  *  process id, and a sequence number. (This is also the behavior if
 323  *  <code>Kernel::rand</code> is called without previously calling
 324  *  <code>srand</code>, but without the sequence.) By setting the seed
 325  *  to a known value, scripts can be made deterministic during testing.
 326  *  The previous seed value is returned. Also see <code>Kernel::rand</code>.
 327  */
 328
 329 static VALUE
 330 rb_f_srand(int argc, VALUE *argv, VALUE obj)
 331 {
 332     VALUE seed, old;
 333
 334     rb_secure(4);
 335     if (argc == 0) {
 336         seed = random_seed();
 337     }
 338     else {
 339         rb_scan_args(argc, argv, "01", &seed);
 340     }
 341     old = rand_init(seed);
 342
 343     return old;
 344 }
 345
 346 static unsigned long
 347 make_mask(unsigned long x)
 348 {
 349     x = x | x >> 1;
 350     x = x | x >> 2;
 351     x = x | x >> 4;
 352     x = x | x >> 8;
 353     x = x | x >> 16;
 354 #if 4 < SIZEOF_LONG
 355     x = x | x >> 32;
 356 #endif
 357     return x;
 358 }
 359
 360 static unsigned long
 361 limited_rand(unsigned long limit)
 362 {
 363     unsigned long mask = make_mask(limit);
 364     int i;
 365     unsigned long val;
 366
 367   retry:
 368     val = 0;
 369     for (i = SIZEOF_LONG/4-1; 0 <= i; i--) {
 370         if (mask >> (i * 32)) {
 371             val |= genrand_int32() << (i * 32);
 372             val &= mask;
 373             if (limit < val)
 374                 goto retry;
 375         }
 376     }
 377     return val;
 378 }
 379
 380 static VALUE
 381 limited_big_rand(struct RBignum *limit)
 382 {
 383     unsigned long mask, lim, rnd;
 384     struct RBignum *val;
 385     int i, len, boundary;
 386
 387     len = (RBIGNUM_LEN(limit) * SIZEOF_BDIGITS + 3) / 4;
 388     val = (struct RBignum *)rb_big_clone((VALUE)limit);
 389     RBIGNUM_SET_SIGN(val, 1);
 390 #if SIZEOF_BDIGITS == 2
 391 # define BIG_GET32(big,i) \
 392     (RBIGNUM_DIGITS(big)[(i)*2] | \
 393      ((i)*2+1 < RBIGNUM_LEN(big) ? \
 394       (RBIGNUM_DIGITS(big)[(i)*2+1] << 16) : \
 395       0))
 396 # define BIG_SET32(big,i,d) \
 397     ((RBIGNUM_DIGITS(big)[(i)*2] = (d) & 0xffff), \
 398      ((i)*2+1 < RBIGNUM_LEN(big) ? \
 399       (RBIGNUM_DIGITS(big)[(i)*2+1] = (d) >> 16) : \
 400       0))
 401 #else
 402     /* SIZEOF_BDIGITS == 4 */
 403 # define BIG_GET32(big,i) (RBIGNUM_DIGITS(big)[i])
 404 # define BIG_SET32(big,i,d) (RBIGNUM_DIGITS(big)[i] = (d))
 405 #endif
 406   retry:
 407     mask = 0;
 408     boundary = 1;
 409     for (i = len-1; 0 <= i; i--) {
 410         lim = BIG_GET32(limit, i);
 411         mask = mask ? 0xffffffff : make_mask(lim);
 412         if (mask) {
 413             rnd = genrand_int32() & mask;
 414             if (boundary) {
 415                 if (lim < rnd)
 416                     goto retry;
 417                 if (rnd < lim)
 418                     boundary = 0;
 419             }
 420         }
 421         else {
 422             rnd = 0;
 423         }
 424         BIG_SET32(val, i, rnd);
 425     }
 426     return rb_big_norm((VALUE)val);
 427 }
 428
 429 /*
 430  *  call-seq:
 431  *     rand(max=0)    => number
 432  *
 433  *  Converts <i>max</i> to an integer using max1 =
 434  *  max<code>.to_i.abs</code>. If the result is zero, returns a
 435  *  pseudorandom floating point number greater than or equal to 0.0 and
 436  *  less than 1.0. Otherwise, returns a pseudorandom integer greater
 437  *  than or equal to zero and less than max1. <code>Kernel::srand</code>
 438  *  may be used to ensure repeatable sequences of random numbers between
 439  *  different runs of the program. Ruby currently uses a modified
 440  *  Mersenne Twister with a period of 2**19937-1.
 441  *
 442  *     srand 1234                 #=> 0
 443  *     [ rand,  rand ]            #=> [0.191519450163469, 0.49766366626136]
 444  *     [ rand(10), rand(1000) ]   #=> [6, 817]
 445  *     srand 1234                 #=> 1234
 446  *     [ rand,  rand ]            #=> [0.191519450163469, 0.49766366626136]
 447  */
 448
 449 static VALUE
 450 rb_f_rand(int argc, VALUE *argv, VALUE obj)
 451 {
 452     VALUE vmax;
 453     long val, max;
 454
 455     rb_scan_args(argc, argv, "01", &vmax);
 456     switch (TYPE(vmax)) {
 457       case T_FLOAT:
 458         if (RFLOAT_VALUE(vmax) <= LONG_MAX && RFLOAT_VALUE(vmax) >= LONG_MIN) {
 459             max = (long)RFLOAT_VALUE(vmax);
 460             break;
 461         }
 462         if (RFLOAT_VALUE(vmax) < 0)
 463             vmax = rb_dbl2big(-RFLOAT_VALUE(vmax));
 464         else
 465             vmax = rb_dbl2big(RFLOAT_VALUE(vmax));
 466         /* fall through */
 467       case T_BIGNUM:
 468       bignum:
 469         {
 470             struct RBignum *limit = (struct RBignum *)vmax;
 471             if (!RBIGNUM_SIGN(limit)) {
 472                 limit = (struct RBignum *)rb_big_clone(vmax);
 473                 RBIGNUM_SET_SIGN(limit, 1);
 474             }
 475             limit = (struct RBignum *)rb_big_minus((VALUE)limit, INT2FIX(1));
 476             if (FIXNUM_P((VALUE)limit)) {
 477                 if (FIX2LONG((VALUE)limit) == -1)
 478                     return DOUBLE2NUM(genrand_real());
 479                 return LONG2NUM(limited_rand(FIX2LONG((VALUE)limit)));
 480             }
 481             return limited_big_rand(limit);
 482         }
 483       case T_NIL:
 484         max = 0;
 485         break;
 486       default:
 487         vmax = rb_Integer(vmax);
 488         if (TYPE(vmax) == T_BIGNUM) goto bignum;
 489       case T_FIXNUM:
 490         max = FIX2LONG(vmax);
 491         break;
 492     }
 493
 494     if (max == 0) {
 495         return DOUBLE2NUM(genrand_real());
 496     }
 497     if (max < 0) max = -max;
 498     val = limited_rand(max-1);
 499     return LONG2NUM(val);
 500 }
 501
 502 void
 503 Init_Random(void)
 504 {
 505     rand_init(random_seed());
 506     rb_define_global_function("srand", rb_f_srand, -1);
 507     rb_define_global_function("rand", rb_f_rand, -1);
 508     rb_global_variable(&saved_seed);
 509 }