• Main Page
  • Modules
  • Data Structures
  • Files
  • File List
  • Globals

random.c

Go to the documentation of this file.
00001 /**********************************************************************
00002 
00003   random.c -
00004 
00005   $Author: nobu $
00006   created at: Fri Dec 24 16:39:21 JST 1993
00007 
00008   Copyright (C) 1993-2007 Yukihiro Matsumoto
00009 
00010 **********************************************************************/
00011 
00012 /*
00013 This is based on trimmed version of MT19937.  To get the original version,
00014 contact <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html>.
00015 
00016 The original copyright notice follows.
00017 
00018    A C-program for MT19937, with initialization improved 2002/2/10.
00019    Coded by Takuji Nishimura and Makoto Matsumoto.
00020    This is a faster version by taking Shawn Cokus's optimization,
00021    Matthe Bellew's simplification, Isaku Wada's real version.
00022 
00023    Before using, initialize the state by using init_genrand(mt, seed)
00024    or init_by_array(mt, init_key, key_length).
00025 
00026    Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
00027    All rights reserved.
00028 
00029    Redistribution and use in source and binary forms, with or without
00030    modification, are permitted provided that the following conditions
00031    are met:
00032 
00033      1. Redistributions of source code must retain the above copyright
00034         notice, this list of conditions and the following disclaimer.
00035 
00036      2. Redistributions in binary form must reproduce the above copyright
00037         notice, this list of conditions and the following disclaimer in the
00038         documentation and/or other materials provided with the distribution.
00039 
00040      3. The names of its contributors may not be used to endorse or promote
00041         products derived from this software without specific prior written
00042         permission.
00043 
00044    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
00045    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
00046    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
00047    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
00048    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00049    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00050    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00051    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00052    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00053    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00054    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00055 
00056 
00057    Any feedback is very welcome.
00058    http://www.math.keio.ac.jp/matumoto/emt.html
00059    email: matumoto@math.keio.ac.jp
00060 */
00061 
00062 #include "ruby/ruby.h"
00063 
00064 #include <limits.h>
00065 #ifdef HAVE_UNISTD_H
00066 #include <unistd.h>
00067 #endif
00068 #include <time.h>
00069 #include <sys/types.h>
00070 #include <sys/stat.h>
00071 #ifdef HAVE_FCNTL_H
00072 #include <fcntl.h>
00073 #endif
00074 #include <math.h>
00075 #include <errno.h>
00076 
00077 #ifdef _WIN32
00078 # if !defined(_WIN32_WINNT) || _WIN32_WINNT < 0x0400
00079 #  undef _WIN32_WINNT
00080 #  define _WIN32_WINNT 0x400
00081 #  undef __WINCRYPT_H__
00082 # endif
00083 #include <wincrypt.h>
00084 #endif
00085 
00086 typedef int int_must_be_32bit_at_least[sizeof(int) * CHAR_BIT < 32 ? -1 : 1];
00087 
00088 /* Period parameters */
00089 #define N 624
00090 #define M 397
00091 #define MATRIX_A 0x9908b0dfU    /* constant vector a */
00092 #define UMASK 0x80000000U       /* most significant w-r bits */
00093 #define LMASK 0x7fffffffU       /* least significant r bits */
00094 #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) )
00095 #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1U ? MATRIX_A : 0U))
00096 
00097 enum {MT_MAX_STATE = N};
00098 
00099 struct MT {
00100     /* assume int is enough to store 32bits */
00101     unsigned int state[N]; /* the array for the state vector  */
00102     unsigned int *next;
00103     int left;
00104 };
00105 
00106 #define genrand_initialized(mt) ((mt)->next != 0)
00107 #define uninit_genrand(mt) ((mt)->next = 0)
00108 
00109 /* initializes state[N] with a seed */
00110 static void
00111 init_genrand(struct MT *mt, unsigned int s)
00112 {
00113     int j;
00114     mt->state[0] = s & 0xffffffffU;
00115     for (j=1; j<N; j++) {
00116         mt->state[j] = (1812433253U * (mt->state[j-1] ^ (mt->state[j-1] >> 30)) + j);
00117         /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
00118         /* In the previous versions, MSBs of the seed affect   */
00119         /* only MSBs of the array state[].                     */
00120         /* 2002/01/09 modified by Makoto Matsumoto             */
00121         mt->state[j] &= 0xffffffff;  /* for >32 bit machines */
00122     }
00123     mt->left = 1;
00124     mt->next = mt->state + N;
00125 }
00126 
00127 /* initialize by an array with array-length */
00128 /* init_key is the array for initializing keys */
00129 /* key_length is its length */
00130 /* slight change for C++, 2004/2/26 */
00131 static void
00132 init_by_array(struct MT *mt, unsigned int init_key[], int key_length)
00133 {
00134     int i, j, k;
00135     init_genrand(mt, 19650218U);
00136     i=1; j=0;
00137     k = (N>key_length ? N : key_length);
00138     for (; k; k--) {
00139         mt->state[i] = (mt->state[i] ^ ((mt->state[i-1] ^ (mt->state[i-1] >> 30)) * 1664525U))
00140           + init_key[j] + j; /* non linear */
00141         mt->state[i] &= 0xffffffffU; /* for WORDSIZE > 32 machines */
00142         i++; j++;
00143         if (i>=N) { mt->state[0] = mt->state[N-1]; i=1; }
00144         if (j>=key_length) j=0;
00145     }
00146     for (k=N-1; k; k--) {
00147         mt->state[i] = (mt->state[i] ^ ((mt->state[i-1] ^ (mt->state[i-1] >> 30)) * 1566083941U))
00148           - i; /* non linear */
00149         mt->state[i] &= 0xffffffffU; /* for WORDSIZE > 32 machines */
00150         i++;
00151         if (i>=N) { mt->state[0] = mt->state[N-1]; i=1; }
00152     }
00153 
00154     mt->state[0] = 0x80000000U; /* MSB is 1; assuring non-zero initial array */
00155 }
00156 
00157 static void
00158 next_state(struct MT *mt)
00159 {
00160     unsigned int *p = mt->state;
00161     int j;
00162 
00163     mt->left = N;
00164     mt->next = mt->state;
00165 
00166     for (j=N-M+1; --j; p++)
00167         *p = p[M] ^ TWIST(p[0], p[1]);
00168 
00169     for (j=M; --j; p++)
00170         *p = p[M-N] ^ TWIST(p[0], p[1]);
00171 
00172     *p = p[M-N] ^ TWIST(p[0], mt->state[0]);
00173 }
00174 
00175 /* generates a random number on [0,0xffffffff]-interval */
00176 static unsigned int
00177 genrand_int32(struct MT *mt)
00178 {
00179     /* mt must be initialized */
00180     unsigned int y;
00181 
00182     if (--mt->left <= 0) next_state(mt);
00183     y = *mt->next++;
00184 
00185     /* Tempering */
00186     y ^= (y >> 11);
00187     y ^= (y << 7) & 0x9d2c5680;
00188     y ^= (y << 15) & 0xefc60000;
00189     y ^= (y >> 18);
00190 
00191     return y;
00192 }
00193 
00194 /* generates a random number on [0,1) with 53-bit resolution*/
00195 static double
00196 genrand_real(struct MT *mt)
00197 {
00198     /* mt must be initialized */
00199     unsigned int a = genrand_int32(mt)>>5, b = genrand_int32(mt)>>6;
00200     return(a*67108864.0+b)*(1.0/9007199254740992.0);
00201 }
00202 
00203 /* generates a random number on [0,1] with 53-bit resolution*/
00204 static double int_pair_to_real_inclusive(unsigned int a, unsigned int b);
00205 static double
00206 genrand_real2(struct MT *mt)
00207 {
00208     /* mt must be initialized */
00209     unsigned int a = genrand_int32(mt), b = genrand_int32(mt);
00210     return int_pair_to_real_inclusive(a, b);
00211 }
00212 
00213 /* These real versions are due to Isaku Wada, 2002/01/09 added */
00214 
00215 #undef N
00216 #undef M
00217 
00218 /* These real versions are due to Isaku Wada, 2002/01/09 added */
00219 
00220 typedef struct {
00221     VALUE seed;
00222     struct MT mt;
00223 } rb_random_t;
00224 
00225 #define DEFAULT_SEED_CNT 4
00226 
00227 static rb_random_t default_rand;
00228 
00229 static VALUE rand_init(struct MT *mt, VALUE vseed);
00230 static VALUE random_seed(void);
00231 
00232 static struct MT *
00233 default_mt(void)
00234 {
00235     rb_random_t *r = &default_rand;
00236     struct MT *mt = &r->mt;
00237     if (!genrand_initialized(mt)) {
00238         r->seed = rand_init(mt, random_seed());
00239     }
00240     return mt;
00241 }
00242 
00243 unsigned int
00244 rb_genrand_int32(void)
00245 {
00246     struct MT *mt = default_mt();
00247     return genrand_int32(mt);
00248 }
00249 
00250 double
00251 rb_genrand_real(void)
00252 {
00253     struct MT *mt = default_mt();
00254     return genrand_real(mt);
00255 }
00256 
00257 #define BDIGITS(x) (RBIGNUM_DIGITS(x))
00258 #define BITSPERDIG (SIZEOF_BDIGITS*CHAR_BIT)
00259 #define BIGRAD ((BDIGIT_DBL)1 << BITSPERDIG)
00260 #define DIGSPERINT (SIZEOF_INT/SIZEOF_BDIGITS)
00261 #define BIGUP(x) ((BDIGIT_DBL)(x) << BITSPERDIG)
00262 #define BIGDN(x) RSHIFT(x,BITSPERDIG)
00263 #define BIGLO(x) ((BDIGIT)((x) & (BIGRAD-1)))
00264 #define BDIGMAX ((BDIGIT)-1)
00265 
00266 #define roomof(n, m) (int)(((n)+(m)-1) / (m))
00267 #define numberof(array) (int)(sizeof(array) / sizeof((array)[0]))
00268 #define SIZEOF_INT32 (31/CHAR_BIT + 1)
00269 
00270 static double
00271 int_pair_to_real_inclusive(unsigned int a, unsigned int b)
00272 {
00273     VALUE x = rb_big_new(roomof(64, BITSPERDIG), 1);
00274     VALUE m = rb_big_new(roomof(53, BITSPERDIG), 1);
00275     BDIGIT *xd = BDIGITS(x);
00276     int i = 0;
00277     double r;
00278 
00279     xd[i++] = (BDIGIT)b;
00280 #if BITSPERDIG < 32
00281     xd[i++] = (BDIGIT)(b >> BITSPERDIG);
00282 #endif
00283     xd[i++] = (BDIGIT)a;
00284 #if BITSPERDIG < 32
00285     xd[i++] = (BDIGIT)(a >> BITSPERDIG);
00286 #endif
00287     xd = BDIGITS(m);
00288 #if BITSPERDIG < 53
00289     MEMZERO(xd, BDIGIT, roomof(53, BITSPERDIG) - 1);
00290 #endif
00291     xd[53 / BITSPERDIG] = 1 << 53 % BITSPERDIG;
00292     xd[0] |= 1;
00293     x = rb_big_mul(x, m);
00294     if (FIXNUM_P(x)) {
00295 #if CHAR_BIT * SIZEOF_LONG > 64
00296         r = (double)(FIX2ULONG(x) >> 64);
00297 #else
00298         return 0.0;
00299 #endif
00300     }
00301     else {
00302 #if 64 % BITSPERDIG == 0
00303         long len = RBIGNUM_LEN(x);
00304         xd = BDIGITS(x);
00305         MEMMOVE(xd, xd + 64 / BITSPERDIG, BDIGIT, len - 64 / BITSPERDIG);
00306         MEMZERO(xd + len - 64 / BITSPERDIG, BDIGIT, 64 / BITSPERDIG);
00307         r = rb_big2dbl(x);
00308 #else
00309         x = rb_big_rshift(x, INT2FIX(64));
00310         if (FIXNUM_P(x)) {
00311             r = (double)FIX2ULONG(x);
00312         }
00313         else {
00314             r = rb_big2dbl(x);
00315         }
00316 #endif
00317     }
00318     return ldexp(r, -53);
00319 }
00320 
00321 VALUE rb_cRandom;
00322 #define id_minus '-'
00323 #define id_plus  '+'
00324 
00325 /* :nodoc: */
00326 static void
00327 random_mark(void *ptr)
00328 {
00329     rb_gc_mark(((rb_random_t *)ptr)->seed);
00330 }
00331 
00332 #define random_free RUBY_TYPED_DEFAULT_FREE
00333 
00334 static size_t
00335 random_memsize(const void *ptr)
00336 {
00337     return ptr ? sizeof(rb_random_t) : 0;
00338 }
00339 
00340 static const rb_data_type_t random_data_type = {
00341     "random",
00342     random_mark,
00343     random_free,
00344     random_memsize,
00345 };
00346 
00347 static rb_random_t *
00348 get_rnd(VALUE obj)
00349 {
00350     rb_random_t *ptr;
00351     TypedData_Get_Struct(obj, rb_random_t, &random_data_type, ptr);
00352     return ptr;
00353 }
00354 
00355 /* :nodoc: */
00356 static VALUE
00357 random_alloc(VALUE klass)
00358 {
00359     rb_random_t *rnd;
00360     VALUE obj = TypedData_Make_Struct(klass, rb_random_t, &random_data_type, rnd);
00361     rnd->seed = INT2FIX(0);
00362     return obj;
00363 }
00364 
00365 static VALUE
00366 rand_init(struct MT *mt, VALUE vseed)
00367 {
00368     volatile VALUE seed;
00369     long blen = 0;
00370     long fixnum_seed;
00371     int i, j, len;
00372     unsigned int buf0[SIZEOF_LONG / SIZEOF_INT32 * 4], *buf = buf0;
00373 
00374     seed = rb_to_int(vseed);
00375     switch (TYPE(seed)) {
00376       case T_FIXNUM:
00377         len = 1;
00378         fixnum_seed = FIX2LONG(seed);
00379         if (fixnum_seed < 0)
00380             fixnum_seed = -fixnum_seed;
00381         buf[0] = (unsigned int)(fixnum_seed & 0xffffffff);
00382 #if SIZEOF_LONG > SIZEOF_INT32
00383         if ((long)(int)fixnum_seed != fixnum_seed) {
00384             if ((buf[1] = (unsigned int)(fixnum_seed >> 32)) != 0) ++len;
00385         }
00386 #endif
00387         break;
00388       case T_BIGNUM:
00389         blen = RBIGNUM_LEN(seed);
00390         if (blen == 0) {
00391             len = 1;
00392         }
00393         else {
00394             if (blen > MT_MAX_STATE * SIZEOF_INT32 / SIZEOF_BDIGITS)
00395                 blen = (len = MT_MAX_STATE) * SIZEOF_INT32 / SIZEOF_BDIGITS;
00396             len = roomof((int)blen * SIZEOF_BDIGITS, SIZEOF_INT32);
00397         }
00398         /* allocate ints for init_by_array */
00399         if (len > numberof(buf0)) buf = ALLOC_N(unsigned int, len);
00400         memset(buf, 0, len * sizeof(*buf));
00401         len = 0;
00402         for (i = (int)(blen-1); 0 <= i; i--) {
00403             j = i * SIZEOF_BDIGITS / SIZEOF_INT32;
00404 #if SIZEOF_BDIGITS < SIZEOF_INT32
00405             buf[j] <<= BITSPERDIG;
00406 #endif
00407             buf[j] |= RBIGNUM_DIGITS(seed)[i];
00408             if (!len && buf[j]) len = j;
00409         }
00410         ++len;
00411         break;
00412       default:
00413         rb_raise(rb_eTypeError, "failed to convert %s into Integer",
00414                  rb_obj_classname(vseed));
00415     }
00416     if (len <= 1) {
00417         init_genrand(mt, buf[0]);
00418     }
00419     else {
00420         if (buf[len-1] == 1) /* remove leading-zero-guard */
00421             len--;
00422         init_by_array(mt, buf, len);
00423     }
00424     if (buf != buf0) xfree(buf);
00425     return seed;
00426 }
00427 
00428 /*
00429  * call-seq: Random.new([seed]) -> prng
00430  *
00431  * Creates new Mersenne Twister based pseudorandom number generator with
00432  * seed.  When the argument seed is omitted, the generator is initialized
00433  * with Random.new_seed.
00434  *
00435  * The argument seed is used to ensure repeatable sequences of random numbers
00436  * between different runs of the program.
00437  *
00438  *     prng = Random.new(1234)
00439  *     [ prng.rand, prng.rand ]   #=> [0.191519450378892, 0.622108771039832]
00440  *     [ prng.integer(10), prng.integer(1000) ]  #=> [4, 664]
00441  *     prng = Random.new(1234)
00442  *     [ prng.rand, prng.rand ]   #=> [0.191519450378892, 0.622108771039832]
00443  */
00444 static VALUE
00445 random_init(int argc, VALUE *argv, VALUE obj)
00446 {
00447     VALUE vseed;
00448     rb_random_t *rnd = get_rnd(obj);
00449 
00450     if (argc == 0) {
00451         vseed = random_seed();
00452     }
00453     else {
00454         rb_scan_args(argc, argv, "01", &vseed);
00455     }
00456     rnd->seed = rand_init(&rnd->mt, vseed);
00457     return obj;
00458 }
00459 
00460 #define DEFAULT_SEED_LEN (DEFAULT_SEED_CNT * sizeof(int))
00461 
00462 #if defined(S_ISCHR) && !defined(DOSISH)
00463 # define USE_DEV_URANDOM 1
00464 #else
00465 # define USE_DEV_URANDOM 0
00466 #endif
00467 
00468 static void
00469 fill_random_seed(unsigned int seed[DEFAULT_SEED_CNT])
00470 {
00471     static int n = 0;
00472     struct timeval tv;
00473 #if USE_DEV_URANDOM
00474     int fd;
00475     struct stat statbuf;
00476 #elif defined(_WIN32)
00477     HCRYPTPROV prov;
00478 #endif
00479 
00480     memset(seed, 0, DEFAULT_SEED_LEN);
00481 
00482 #if USE_DEV_URANDOM
00483     if ((fd = open("/dev/urandom", O_RDONLY
00484 #ifdef O_NONBLOCK
00485             |O_NONBLOCK
00486 #endif
00487 #ifdef O_NOCTTY
00488             |O_NOCTTY
00489 #endif
00490             )) >= 0) {
00491         if (fstat(fd, &statbuf) == 0 && S_ISCHR(statbuf.st_mode)) {
00492             (void)read(fd, seed, DEFAULT_SEED_LEN);
00493         }
00494         close(fd);
00495     }
00496 #elif defined(_WIN32)
00497     if (CryptAcquireContext(&prov, NULL, NULL, PROV_RSA_FULL, CRYPT_VERIFYCONTEXT)) {
00498         CryptGenRandom(prov, DEFAULT_SEED_LEN, (void *)seed);
00499         CryptReleaseContext(prov, 0);
00500     }
00501 #endif
00502 
00503     gettimeofday(&tv, 0);
00504     seed[0] ^= tv.tv_usec;
00505     seed[1] ^= (unsigned int)tv.tv_sec;
00506 #if SIZEOF_TIME_T > SIZEOF_INT
00507     seed[0] ^= (unsigned int)((time_t)tv.tv_sec >> SIZEOF_INT * CHAR_BIT);
00508 #endif
00509     seed[2] ^= getpid() ^ (n++ << 16);
00510     seed[3] ^= (unsigned int)(VALUE)&seed;
00511 #if SIZEOF_VOIDP > SIZEOF_INT
00512     seed[2] ^= (unsigned int)((VALUE)&seed >> SIZEOF_INT * CHAR_BIT);
00513 #endif
00514 }
00515 
00516 static VALUE
00517 make_seed_value(const void *ptr)
00518 {
00519     const long len = DEFAULT_SEED_LEN/SIZEOF_BDIGITS;
00520     BDIGIT *digits;
00521     NEWOBJ(big, struct RBignum);
00522     OBJSETUP(big, rb_cBignum, T_BIGNUM);
00523 
00524     RBIGNUM_SET_SIGN(big, 1);
00525     rb_big_resize((VALUE)big, len + 1);
00526     digits = RBIGNUM_DIGITS(big);
00527 
00528     MEMCPY(digits, ptr, char, DEFAULT_SEED_LEN);
00529 
00530     /* set leading-zero-guard if need. */
00531     digits[len] =
00532 #if SIZEOF_INT32 / SIZEOF_BDIGITS > 1
00533         digits[len-2] <= 1 && digits[len-1] == 0
00534 #else
00535         digits[len-1] <= 1
00536 #endif
00537         ? 1 : 0;
00538 
00539     return rb_big_norm((VALUE)big);
00540 }
00541 
00542 /*
00543  * call-seq: Random.new_seed -> integer
00544  *
00545  * Returns arbitrary value for seed.
00546  */
00547 static VALUE
00548 random_seed(void)
00549 {
00550     unsigned int buf[DEFAULT_SEED_CNT];
00551     fill_random_seed(buf);
00552     return make_seed_value(buf);
00553 }
00554 
00555 /*
00556  * call-seq: prng.seed -> integer
00557  *
00558  * Returns the seed of the generator.
00559  */
00560 static VALUE
00561 random_get_seed(VALUE obj)
00562 {
00563     return get_rnd(obj)->seed;
00564 }
00565 
00566 /* :nodoc: */
00567 static VALUE
00568 random_copy(VALUE obj, VALUE orig)
00569 {
00570     rb_random_t *rnd1 = get_rnd(obj);
00571     rb_random_t *rnd2 = get_rnd(orig);
00572     struct MT *mt = &rnd1->mt;
00573 
00574     *rnd1 = *rnd2;
00575     mt->next = mt->state + numberof(mt->state) - mt->left + 1;
00576     return obj;
00577 }
00578 
00579 static VALUE
00580 mt_state(const struct MT *mt)
00581 {
00582     VALUE bigo = rb_big_new(sizeof(mt->state) / sizeof(BDIGIT), 1);
00583     BDIGIT *d = RBIGNUM_DIGITS(bigo);
00584     int i;
00585 
00586     for (i = 0; i < numberof(mt->state); ++i) {
00587         unsigned int x = mt->state[i];
00588 #if SIZEOF_BDIGITS < SIZEOF_INT32
00589         int j;
00590         for (j = 0; j < SIZEOF_INT32 / SIZEOF_BDIGITS; ++j) {
00591             *d++ = BIGLO(x);
00592             x = BIGDN(x);
00593         }
00594 #else
00595         *d++ = (BDIGIT)x;
00596 #endif
00597     }
00598     return rb_big_norm(bigo);
00599 }
00600 
00601 /* :nodoc: */
00602 static VALUE
00603 random_state(VALUE obj)
00604 {
00605     rb_random_t *rnd = get_rnd(obj);
00606     return mt_state(&rnd->mt);
00607 }
00608 
00609 /* :nodoc: */
00610 static VALUE
00611 random_s_state(VALUE klass)
00612 {
00613     return mt_state(&default_rand.mt);
00614 }
00615 
00616 /* :nodoc: */
00617 static VALUE
00618 random_left(VALUE obj)
00619 {
00620     rb_random_t *rnd = get_rnd(obj);
00621     return INT2FIX(rnd->mt.left);
00622 }
00623 
00624 /* :nodoc: */
00625 static VALUE
00626 random_s_left(VALUE klass)
00627 {
00628     return INT2FIX(default_rand.mt.left);
00629 }
00630 
00631 /* :nodoc: */
00632 static VALUE
00633 random_dump(VALUE obj)
00634 {
00635     rb_random_t *rnd = get_rnd(obj);
00636     VALUE dump = rb_ary_new2(3);
00637 
00638     rb_ary_push(dump, mt_state(&rnd->mt));
00639     rb_ary_push(dump, INT2FIX(rnd->mt.left));
00640     rb_ary_push(dump, rnd->seed);
00641 
00642     return dump;
00643 }
00644 
00645 /* :nodoc: */
00646 static VALUE
00647 random_load(VALUE obj, VALUE dump)
00648 {
00649     rb_random_t *rnd = get_rnd(obj);
00650     struct MT *mt = &rnd->mt;
00651     VALUE state, left = INT2FIX(1), seed = INT2FIX(0);
00652     VALUE *ary;
00653     unsigned long x;
00654 
00655     Check_Type(dump, T_ARRAY);
00656     ary = RARRAY_PTR(dump);
00657     switch (RARRAY_LEN(dump)) {
00658       case 3:
00659         seed = ary[2];
00660       case 2:
00661         left = ary[1];
00662       case 1:
00663         state = ary[0];
00664         break;
00665       default:
00666         rb_raise(rb_eArgError, "wrong dump data");
00667     }
00668     memset(mt->state, 0, sizeof(mt->state));
00669     if (FIXNUM_P(state)) {
00670         x = FIX2ULONG(state);
00671         mt->state[0] = (unsigned int)x;
00672 #if SIZEOF_LONG / SIZEOF_INT >= 2
00673         mt->state[1] = (unsigned int)(x >> BITSPERDIG);
00674 #endif
00675 #if SIZEOF_LONG / SIZEOF_INT >= 3
00676         mt->state[2] = (unsigned int)(x >> 2 * BITSPERDIG);
00677 #endif
00678 #if SIZEOF_LONG / SIZEOF_INT >= 4
00679         mt->state[3] = (unsigned int)(x >> 3 * BITSPERDIG);
00680 #endif
00681     }
00682     else {
00683         BDIGIT *d;
00684         long len;
00685         Check_Type(state, T_BIGNUM);
00686         len = RBIGNUM_LEN(state);
00687         if (len > roomof(sizeof(mt->state), SIZEOF_BDIGITS)) {
00688             len = roomof(sizeof(mt->state), SIZEOF_BDIGITS);
00689         }
00690 #if SIZEOF_BDIGITS < SIZEOF_INT
00691         else if (len % DIGSPERINT) {
00692             d = RBIGNUM_DIGITS(state) + len;
00693 # if DIGSPERINT == 2
00694             --len;
00695             x = *--d;
00696 # else
00697             x = 0;
00698             do {
00699                 x = (x << BITSPERDIG) | *--d;
00700             } while (--len % DIGSPERINT);
00701 # endif
00702             mt->state[len / DIGSPERINT] = (unsigned int)x;
00703         }
00704 #endif
00705         if (len > 0) {
00706             d = BDIGITS(state) + len;
00707             do {
00708                 --len;
00709                 x = *--d;
00710 # if DIGSPERINT == 2
00711                 --len;
00712                 x = (x << BITSPERDIG) | *--d;
00713 # elif SIZEOF_BDIGITS < SIZEOF_INT
00714                 do {
00715                     x = (x << BITSPERDIG) | *--d;
00716                 } while (--len % DIGSPERINT);
00717 # endif
00718                 mt->state[len / DIGSPERINT] = (unsigned int)x;
00719             } while (len > 0);
00720         }
00721     }
00722     x = NUM2ULONG(left);
00723     if (x > numberof(mt->state)) {
00724         rb_raise(rb_eArgError, "wrong value");
00725     }
00726     mt->left = (unsigned int)x;
00727     mt->next = mt->state + numberof(mt->state) - x + 1;
00728     rnd->seed = rb_to_int(seed);
00729 
00730     return obj;
00731 }
00732 
00733 /*
00734  *  call-seq:
00735  *     srand(number=0)    -> old_seed
00736  *
00737  *  Seeds the pseudorandom number generator to the value of
00738  *  <i>number</i>. If <i>number</i> is omitted
00739  *  or zero, seeds the generator using a combination of the time, the
00740  *  process id, and a sequence number. (This is also the behavior if
00741  *  <code>Kernel::rand</code> is called without previously calling
00742  *  <code>srand</code>, but without the sequence.) By setting the seed
00743  *  to a known value, scripts can be made deterministic during testing.
00744  *  The previous seed value is returned. Also see <code>Kernel::rand</code>.
00745  */
00746 
00747 static VALUE
00748 rb_f_srand(int argc, VALUE *argv, VALUE obj)
00749 {
00750     VALUE seed, old;
00751     rb_random_t *r = &default_rand;
00752 
00753     rb_secure(4);
00754     if (argc == 0) {
00755         seed = random_seed();
00756     }
00757     else {
00758         rb_scan_args(argc, argv, "01", &seed);
00759     }
00760     old = r->seed;
00761     r->seed = rand_init(&r->mt, seed);
00762 
00763     return old;
00764 }
00765 
00766 static unsigned long
00767 make_mask(unsigned long x)
00768 {
00769     x = x | x >> 1;
00770     x = x | x >> 2;
00771     x = x | x >> 4;
00772     x = x | x >> 8;
00773     x = x | x >> 16;
00774 #if 4 < SIZEOF_LONG
00775     x = x | x >> 32;
00776 #endif
00777     return x;
00778 }
00779 
00780 static unsigned long
00781 limited_rand(struct MT *mt, unsigned long limit)
00782 {
00783     /* mt must be initialized */
00784     int i;
00785     unsigned long val, mask;
00786 
00787     if (!limit) return 0;
00788     mask = make_mask(limit);
00789   retry:
00790     val = 0;
00791     for (i = SIZEOF_LONG/SIZEOF_INT32-1; 0 <= i; i--) {
00792         if ((mask >> (i * 32)) & 0xffffffff) {
00793             val |= (unsigned long)genrand_int32(mt) << (i * 32);
00794             val &= mask;
00795             if (limit < val)
00796                 goto retry;
00797         }
00798     }
00799     return val;
00800 }
00801 
00802 static VALUE
00803 limited_big_rand(struct MT *mt, struct RBignum *limit)
00804 {
00805     /* mt must be initialized */
00806     unsigned long mask, lim, rnd;
00807     struct RBignum *val;
00808     long i, len;
00809     int boundary;
00810 
00811     len = (RBIGNUM_LEN(limit) * SIZEOF_BDIGITS + 3) / 4;
00812     val = (struct RBignum *)rb_big_clone((VALUE)limit);
00813     RBIGNUM_SET_SIGN(val, 1);
00814 #if SIZEOF_BDIGITS == 2
00815 # define BIG_GET32(big,i) \
00816     (RBIGNUM_DIGITS(big)[(i)*2] | \
00817      ((i)*2+1 < RBIGNUM_LEN(big) ? \
00818       (RBIGNUM_DIGITS(big)[(i)*2+1] << 16) : \
00819       0))
00820 # define BIG_SET32(big,i,d) \
00821     ((RBIGNUM_DIGITS(big)[(i)*2] = (d) & 0xffff), \
00822      ((i)*2+1 < RBIGNUM_LEN(big) ? \
00823       (RBIGNUM_DIGITS(big)[(i)*2+1] = (d) >> 16) : \
00824       0))
00825 #else
00826     /* SIZEOF_BDIGITS == 4 */
00827 # define BIG_GET32(big,i) (RBIGNUM_DIGITS(big)[i])
00828 # define BIG_SET32(big,i,d) (RBIGNUM_DIGITS(big)[i] = (d))
00829 #endif
00830   retry:
00831     mask = 0;
00832     boundary = 1;
00833     for (i = len-1; 0 <= i; i--) {
00834         lim = BIG_GET32(limit, i);
00835         mask = mask ? 0xffffffff : make_mask(lim);
00836         if (mask) {
00837             rnd = genrand_int32(mt) & mask;
00838             if (boundary) {
00839                 if (lim < rnd)
00840                     goto retry;
00841                 if (rnd < lim)
00842                     boundary = 0;
00843             }
00844         }
00845         else {
00846             rnd = 0;
00847         }
00848         BIG_SET32(val, i, (BDIGIT)rnd);
00849     }
00850     return rb_big_norm((VALUE)val);
00851 }
00852 
00853 unsigned long
00854 rb_rand_internal(unsigned long i)
00855 {
00856     struct MT *mt = default_mt();
00857     return limited_rand(mt, i);
00858 }
00859 
00860 unsigned int
00861 rb_random_int32(VALUE obj)
00862 {
00863     rb_random_t *rnd = get_rnd(obj);
00864     return genrand_int32(&rnd->mt);
00865 }
00866 
00867 double
00868 rb_random_real(VALUE obj)
00869 {
00870     rb_random_t *rnd = get_rnd(obj);
00871     return genrand_real(&rnd->mt);
00872 }
00873 
00874 /*
00875  * call-seq: prng.bytes(size) -> prng
00876  *
00877  * Returns a random binary string.  The argument size specified the length of
00878  * the result string.
00879  */
00880 static VALUE
00881 random_bytes(VALUE obj, VALUE len)
00882 {
00883     return rb_random_bytes(obj, NUM2LONG(rb_to_int(len)));
00884 }
00885 
00886 VALUE
00887 rb_random_bytes(VALUE obj, long n)
00888 {
00889     rb_random_t *rnd = get_rnd(obj);
00890     VALUE bytes = rb_str_new(0, n);
00891     char *ptr = RSTRING_PTR(bytes);
00892     unsigned int r, i;
00893 
00894     for (; n >= SIZEOF_INT32; n -= SIZEOF_INT32) {
00895         r = genrand_int32(&rnd->mt);
00896         i = SIZEOF_INT32;
00897         do {
00898             *ptr++ = (char)r;
00899             r >>= CHAR_BIT;
00900         } while (--i);
00901     }
00902     if (n > 0) {
00903         r = genrand_int32(&rnd->mt);
00904         do {
00905             *ptr++ = (char)r;
00906             r >>= CHAR_BIT;
00907         } while (--n);
00908     }
00909     return bytes;
00910 }
00911 
00912 static VALUE
00913 range_values(VALUE vmax, VALUE *begp, int *exclp)
00914 {
00915     VALUE end, r;
00916 
00917     if (!rb_range_values(vmax, begp, &end, exclp)) return Qfalse;
00918     if (!rb_respond_to(end, id_minus)) return Qfalse;
00919     r = rb_funcall2(end, id_minus, 1, begp);
00920     if (NIL_P(r)) return Qfalse;
00921     return r;
00922 }
00923 
00924 static VALUE
00925 rand_int(struct MT *mt, VALUE vmax, int restrictive)
00926 {
00927     /* mt must be initialized */
00928     long max;
00929     unsigned long r;
00930 
00931     if (FIXNUM_P(vmax)) {
00932         max = FIX2LONG(vmax);
00933         if (!max) return Qnil;
00934         if (max < 0) {
00935             if (restrictive) return Qnil;
00936             max = -max;
00937         }
00938         r = limited_rand(mt, (unsigned long)max - 1);
00939         return ULONG2NUM(r);
00940     }
00941     else {
00942         VALUE ret;
00943         if (rb_bigzero_p(vmax)) return Qnil;
00944         if (!RBIGNUM_SIGN(vmax)) {
00945             if (restrictive) return Qnil;
00946             vmax = rb_big_clone(vmax);
00947             RBIGNUM_SET_SIGN(vmax, 1);
00948         }
00949         vmax = rb_big_minus(vmax, INT2FIX(1));
00950         if (FIXNUM_P(vmax)) {
00951             max = FIX2LONG(vmax);
00952             if (max == -1) return Qnil;
00953             r = limited_rand(mt, max);
00954             return LONG2NUM(r);
00955         }
00956         ret = limited_big_rand(mt, RBIGNUM(vmax));
00957         RB_GC_GUARD(vmax);
00958         return ret;
00959     }
00960 }
00961 
00962 static inline double
00963 float_value(VALUE v)
00964 {
00965     double x = RFLOAT_VALUE(v);
00966     if (isinf(x) || isnan(x)) {
00967         VALUE error = INT2FIX(EDOM);
00968         rb_exc_raise(rb_class_new_instance(1, &error, rb_eSystemCallError));
00969     }
00970     return x;
00971 }
00972 
00973 /*
00974  * call-seq:
00975  *     prng.rand -> float
00976  *     prng.rand(limit) -> number
00977  *
00978  * When the argument is an +Integer+ or a +Bignum+, it returns a
00979  * random integer greater than or equal to zero and less than the
00980  * argument.  Unlike Random.rand, when the argument is a negative
00981  * integer or zero, it raises an ArgumentError.
00982  *
00983  * When the argument is a +Float+, it returns a random floating point
00984  * number between 0.0 and _max_, including 0.0 and excluding _max_.
00985  *
00986  * When the argument _limit_ is a +Range+, it returns a random
00987  * number where range.member?(number) == true.
00988  *     prng.rand(5..9)  #=> one of [5, 6, 7, 8, 9]
00989  *     prng.rand(5...9) #=> one of [5, 6, 7, 8]
00990  *     prng.rand(5.0..9.0) #=> between 5.0 and 9.0, including 9.0
00991  *     prng.rand(5.0...9.0) #=> between 5.0 and 9.0, excluding 9.0
00992  *
00993  * +begin+/+end+ of the range have to have subtract and add methods.
00994  *
00995  * Otherwise, it raises an ArgumentError.
00996  */
00997 static VALUE
00998 random_rand(int argc, VALUE *argv, VALUE obj)
00999 {
01000     rb_random_t *rnd = get_rnd(obj);
01001     VALUE vmax, beg = Qundef, v;
01002     int excl = 0;
01003 
01004     if (argc == 0) {
01005         return rb_float_new(genrand_real(&rnd->mt));
01006     }
01007     else if (argc != 1) {
01008         rb_raise(rb_eArgError, "wrong number of arguments (%d for 0..1)", argc);
01009     }
01010     vmax = argv[0];
01011     if (NIL_P(vmax)) {
01012         v = Qnil;
01013     }
01014     else if (TYPE(vmax) != T_FLOAT && (v = rb_check_to_integer(vmax, "to_int"), !NIL_P(v))) {
01015         v = rand_int(&rnd->mt, vmax = v, 1);
01016     }
01017     else if (v = rb_check_to_float(vmax), !NIL_P(v)) {
01018         double max = float_value(v);
01019         if (max > 0.0)
01020             v = rb_float_new(max * genrand_real(&rnd->mt));
01021         else
01022             v = Qnil;
01023     }
01024     else if ((v = range_values(vmax, &beg, &excl)) != Qfalse) {
01025         vmax = v;
01026         if (TYPE(vmax) != T_FLOAT && (v = rb_check_to_integer(vmax, "to_int"), !NIL_P(v))) {
01027             long max;
01028             vmax = v;
01029             v = Qnil;
01030             if (FIXNUM_P(vmax)) {
01031               fixnum:
01032                 if ((max = FIX2LONG(vmax) - excl) >= 0) {
01033                     unsigned long r = limited_rand(&rnd->mt, (unsigned long)max);
01034                     v = ULONG2NUM(r);
01035                 }
01036             }
01037             else if (BUILTIN_TYPE(vmax) == T_BIGNUM && RBIGNUM_SIGN(vmax) && !rb_bigzero_p(vmax)) {
01038                 vmax = excl ? rb_big_minus(vmax, INT2FIX(1)) : rb_big_norm(vmax);
01039                 if (FIXNUM_P(vmax)) {
01040                     excl = 0;
01041                     goto fixnum;
01042                 }
01043                 v = limited_big_rand(&rnd->mt, RBIGNUM(vmax));
01044             }
01045         }
01046         else if (v = rb_check_to_float(vmax), !NIL_P(v)) {
01047             double max = float_value(v), r;
01048             v = Qnil;
01049             if (max > 0.0) {
01050                 if (excl) {
01051                     r = genrand_real(&rnd->mt);
01052                 }
01053                 else {
01054                     r = genrand_real2(&rnd->mt);
01055                 }
01056                 v = rb_float_new(r * max);
01057             }
01058             else if (max == 0.0 && !excl) {
01059                 v = rb_float_new(0.0);
01060             }
01061         }
01062     }
01063     else {
01064         v = Qnil;
01065         NUM2LONG(vmax);
01066     }
01067     if (NIL_P(v)) {
01068         VALUE mesg = rb_str_new_cstr("invalid argument - ");
01069         rb_str_append(mesg, rb_obj_as_string(argv[0]));
01070         rb_exc_raise(rb_exc_new3(rb_eArgError, mesg));
01071     }
01072     if (beg == Qundef) return v;
01073     if (FIXNUM_P(beg) && FIXNUM_P(v)) {
01074         long x = FIX2LONG(beg) + FIX2LONG(v);
01075         return LONG2NUM(x);
01076     }
01077     switch (TYPE(v)) {
01078       case T_BIGNUM:
01079         return rb_big_plus(v, beg);
01080       case T_FLOAT: {
01081         VALUE f = rb_check_to_float(beg);
01082         if (!NIL_P(f)) {
01083             RFLOAT_VALUE(v) += RFLOAT_VALUE(f);
01084             return v;
01085         }
01086       }
01087       default:
01088         return rb_funcall2(beg, id_plus, 1, &v);
01089     }
01090 }
01091 
01092 /*
01093  * call-seq:
01094  *     prng1 == prng2 -> true or false
01095  *
01096  * Returns true if the generators' states equal.
01097  */
01098 static VALUE
01099 random_equal(VALUE self, VALUE other)
01100 {
01101     rb_random_t *r1, *r2;
01102     if (rb_obj_class(self) != rb_obj_class(other)) return Qfalse;
01103     r1 = get_rnd(self);
01104     r2 = get_rnd(other);
01105     if (!RTEST(rb_funcall2(r1->seed, rb_intern("=="), 1, &r2->seed))) return Qfalse;
01106     if (memcmp(r1->mt.state, r2->mt.state, sizeof(r1->mt.state))) return Qfalse;
01107     if ((r1->mt.next - r1->mt.state) != (r2->mt.next - r2->mt.state)) return Qfalse;
01108     if (r1->mt.left != r2->mt.left) return Qfalse;
01109     return Qtrue;
01110 }
01111 
01112 /*
01113  *  call-seq:
01114  *     rand(max=0)    -> number
01115  *
01116  *  Converts <i>max</i> to an integer using max1 =
01117  *  max<code>.to_i.abs</code>. If _max_ is +nil+ the result is zero, returns a
01118  *  pseudorandom floating point number greater than or equal to 0.0 and
01119  *  less than 1.0. Otherwise, returns a pseudorandom integer greater
01120  *  than or equal to zero and less than max1. <code>Kernel::srand</code>
01121  *  may be used to ensure repeatable sequences of random numbers between
01122  *  different runs of the program. Ruby currently uses a modified
01123  *  Mersenne Twister with a period of 2**19937-1.
01124  *
01125  *     srand 1234                 #=> 0
01126  *     [ rand,  rand ]            #=> [0.191519450163469, 0.49766366626136]
01127  *     [ rand(10), rand(1000) ]   #=> [6, 817]
01128  *     srand 1234                 #=> 1234
01129  *     [ rand,  rand ]            #=> [0.191519450163469, 0.49766366626136]
01130  */
01131 
01132 static VALUE
01133 rb_f_rand(int argc, VALUE *argv, VALUE obj)
01134 {
01135     VALUE vmax, r;
01136     struct MT *mt = default_mt();
01137 
01138     if (argc == 0) goto zero_arg;
01139     rb_scan_args(argc, argv, "01", &vmax);
01140     if (NIL_P(vmax)) goto zero_arg;
01141     vmax = rb_to_int(vmax);
01142     if (vmax == INT2FIX(0) || NIL_P(r = rand_int(mt, vmax, 0))) {
01143       zero_arg:
01144         return DBL2NUM(genrand_real(mt));
01145     }
01146     return r;
01147 }
01148 
01149 static st_index_t hashseed;
01150 
01151 static VALUE
01152 init_randomseed(struct MT *mt, unsigned int initial[DEFAULT_SEED_CNT])
01153 {
01154     VALUE seed;
01155     fill_random_seed(initial);
01156     init_by_array(mt, initial, DEFAULT_SEED_CNT);
01157     seed = make_seed_value(initial);
01158     memset(initial, 0, DEFAULT_SEED_LEN);
01159     return seed;
01160 }
01161 
01162 void
01163 Init_RandomSeed(void)
01164 {
01165     rb_random_t *r = &default_rand;
01166     unsigned int initial[DEFAULT_SEED_CNT];
01167     struct MT *mt = &r->mt;
01168     VALUE seed = init_randomseed(mt, initial);
01169 
01170     hashseed = genrand_int32(mt);
01171 #if SIZEOF_ST_INDEX_T*CHAR_BIT > 4*8
01172     hashseed <<= 32;
01173     hashseed |= genrand_int32(mt);
01174 #endif
01175 #if SIZEOF_ST_INDEX_T*CHAR_BIT > 8*8
01176     hashseed <<= 32;
01177     hashseed |= genrand_int32(mt);
01178 #endif
01179 #if SIZEOF_ST_INDEX_T*CHAR_BIT > 12*8
01180     hashseed <<= 32;
01181     hashseed |= genrand_int32(mt);
01182 #endif
01183 
01184     rb_global_variable(&r->seed);
01185     r->seed = seed;
01186 }
01187 
01188 st_index_t
01189 rb_hash_start(st_index_t h)
01190 {
01191     return st_hash_start(hashseed + h);
01192 }
01193 
01194 static void
01195 Init_RandomSeed2(void)
01196 {
01197     VALUE seed = default_rand.seed;
01198 
01199     if (RB_TYPE_P(seed, T_BIGNUM)) {
01200         RBASIC(seed)->klass = rb_cBignum;
01201     }
01202 }
01203 
01204 void
01205 rb_reset_random_seed(void)
01206 {
01207     rb_random_t *r = &default_rand;
01208     uninit_genrand(&r->mt);
01209     r->seed = INT2FIX(0);
01210 }
01211 
01212 void
01213 Init_Random(void)
01214 {
01215     Init_RandomSeed2();
01216     rb_define_global_function("srand", rb_f_srand, -1);
01217     rb_define_global_function("rand", rb_f_rand, -1);
01218 
01219     rb_cRandom = rb_define_class("Random", rb_cObject);
01220     rb_define_alloc_func(rb_cRandom, random_alloc);
01221     rb_define_method(rb_cRandom, "initialize", random_init, -1);
01222     rb_define_method(rb_cRandom, "rand", random_rand, -1);
01223     rb_define_method(rb_cRandom, "bytes", random_bytes, 1);
01224     rb_define_method(rb_cRandom, "seed", random_get_seed, 0);
01225     rb_define_method(rb_cRandom, "initialize_copy", random_copy, 1);
01226     rb_define_method(rb_cRandom, "marshal_dump", random_dump, 0);
01227     rb_define_method(rb_cRandom, "marshal_load", random_load, 1);
01228     rb_define_private_method(rb_cRandom, "state", random_state, 0);
01229     rb_define_private_method(rb_cRandom, "left", random_left, 0);
01230     rb_define_method(rb_cRandom, "==", random_equal, 1);
01231 
01232     rb_define_singleton_method(rb_cRandom, "srand", rb_f_srand, -1);
01233     rb_define_singleton_method(rb_cRandom, "rand", rb_f_rand, -1);
01234     rb_define_singleton_method(rb_cRandom, "new_seed", random_seed, 0);
01235     rb_define_private_method(CLASS_OF(rb_cRandom), "state", random_s_state, 0);
01236     rb_define_private_method(CLASS_OF(rb_cRandom), "left", random_s_left, 0);
01237 }
01238 

Generated on Wed Sep 8 2010 21:55:11 for Ruby by  doxygen 1.7.1