Ruby  1.9.3p429(2013-05-15revision40747)
random.c
Go to the documentation of this file.
1 /**********************************************************************
2 
3  random.c -
4 
5  $Author: usa $
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(mt, seed)
24  or init_by_array(mt, 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 #include "ruby/ruby.h"
63 
64 #include <limits.h>
65 #ifdef HAVE_UNISTD_H
66 #include <unistd.h>
67 #endif
68 #include <time.h>
69 #include <sys/types.h>
70 #include <sys/stat.h>
71 #ifdef HAVE_FCNTL_H
72 #include <fcntl.h>
73 #endif
74 #include <math.h>
75 #include <errno.h>
76 #if defined(HAVE_SYS_TIME_H)
77 #include <sys/time.h>
78 #endif
79 
80 #ifdef _WIN32
81 # if !defined(_WIN32_WINNT) || _WIN32_WINNT < 0x0400
82 # undef _WIN32_WINNT
83 # define _WIN32_WINNT 0x400
84 # undef __WINCRYPT_H__
85 # endif
86 #include <wincrypt.h>
87 #endif
88 
89 typedef int int_must_be_32bit_at_least[sizeof(int) * CHAR_BIT < 32 ? -1 : 1];
90 
91 /* Period parameters */
92 #define N 624
93 #define M 397
94 #define MATRIX_A 0x9908b0dfU /* constant vector a */
95 #define UMASK 0x80000000U /* most significant w-r bits */
96 #define LMASK 0x7fffffffU /* least significant r bits */
97 #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) )
98 #define TWIST(u,v) ((MIXBITS((u),(v)) >> 1) ^ ((v)&1U ? MATRIX_A : 0U))
99 
100 enum {MT_MAX_STATE = N};
101 
102 struct MT {
103  /* assume int is enough to store 32bits */
104  unsigned int state[N]; /* the array for the state vector */
105  unsigned int *next;
106  int left;
107 };
108 
109 #define genrand_initialized(mt) ((mt)->next != 0)
110 #define uninit_genrand(mt) ((mt)->next = 0)
111 
112 /* initializes state[N] with a seed */
113 static void
114 init_genrand(struct MT *mt, unsigned int s)
115 {
116  int j;
117  mt->state[0] = s & 0xffffffffU;
118  for (j=1; j<N; j++) {
119  mt->state[j] = (1812433253U * (mt->state[j-1] ^ (mt->state[j-1] >> 30)) + j);
120  /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
121  /* In the previous versions, MSBs of the seed affect */
122  /* only MSBs of the array state[]. */
123  /* 2002/01/09 modified by Makoto Matsumoto */
124  mt->state[j] &= 0xffffffff; /* for >32 bit machines */
125  }
126  mt->left = 1;
127  mt->next = mt->state + N;
128 }
129 
130 /* initialize by an array with array-length */
131 /* init_key is the array for initializing keys */
132 /* key_length is its length */
133 /* slight change for C++, 2004/2/26 */
134 static void
135 init_by_array(struct MT *mt, unsigned int init_key[], int key_length)
136 {
137  int i, j, k;
138  init_genrand(mt, 19650218U);
139  i=1; j=0;
140  k = (N>key_length ? N : key_length);
141  for (; k; k--) {
142  mt->state[i] = (mt->state[i] ^ ((mt->state[i-1] ^ (mt->state[i-1] >> 30)) * 1664525U))
143  + init_key[j] + j; /* non linear */
144  mt->state[i] &= 0xffffffffU; /* for WORDSIZE > 32 machines */
145  i++; j++;
146  if (i>=N) { mt->state[0] = mt->state[N-1]; i=1; }
147  if (j>=key_length) j=0;
148  }
149  for (k=N-1; k; k--) {
150  mt->state[i] = (mt->state[i] ^ ((mt->state[i-1] ^ (mt->state[i-1] >> 30)) * 1566083941U))
151  - i; /* non linear */
152  mt->state[i] &= 0xffffffffU; /* for WORDSIZE > 32 machines */
153  i++;
154  if (i>=N) { mt->state[0] = mt->state[N-1]; i=1; }
155  }
156 
157  mt->state[0] = 0x80000000U; /* MSB is 1; assuring non-zero initial array */
158 }
159 
160 static void
161 next_state(struct MT *mt)
162 {
163  unsigned int *p = mt->state;
164  int j;
165 
166  mt->left = N;
167  mt->next = mt->state;
168 
169  for (j=N-M+1; --j; p++)
170  *p = p[M] ^ TWIST(p[0], p[1]);
171 
172  for (j=M; --j; p++)
173  *p = p[M-N] ^ TWIST(p[0], p[1]);
174 
175  *p = p[M-N] ^ TWIST(p[0], mt->state[0]);
176 }
177 
178 /* generates a random number on [0,0xffffffff]-interval */
179 static unsigned int
180 genrand_int32(struct MT *mt)
181 {
182  /* mt must be initialized */
183  unsigned int y;
184 
185  if (--mt->left <= 0) next_state(mt);
186  y = *mt->next++;
187 
188  /* Tempering */
189  y ^= (y >> 11);
190  y ^= (y << 7) & 0x9d2c5680;
191  y ^= (y << 15) & 0xefc60000;
192  y ^= (y >> 18);
193 
194  return y;
195 }
196 
197 /* generates a random number on [0,1) with 53-bit resolution*/
198 static double
199 genrand_real(struct MT *mt)
200 {
201  /* mt must be initialized */
202  unsigned int a = genrand_int32(mt)>>5, b = genrand_int32(mt)>>6;
203  return(a*67108864.0+b)*(1.0/9007199254740992.0);
204 }
205 
206 /* generates a random number on [0,1] with 53-bit resolution*/
207 static double int_pair_to_real_inclusive(unsigned int a, unsigned int b);
208 static double
209 genrand_real2(struct MT *mt)
210 {
211  /* mt must be initialized */
212  unsigned int a = genrand_int32(mt), b = genrand_int32(mt);
213  return int_pair_to_real_inclusive(a, b);
214 }
215 
216 /* These real versions are due to Isaku Wada, 2002/01/09 added */
217 
218 #undef N
219 #undef M
220 
221 /* These real versions are due to Isaku Wada, 2002/01/09 added */
222 
223 typedef struct {
225  struct MT mt;
226 } rb_random_t;
227 
228 #define DEFAULT_SEED_CNT 4
229 
231 
232 static VALUE rand_init(struct MT *mt, VALUE vseed);
233 static VALUE random_seed(void);
234 
235 static rb_random_t *
237 {
238  struct MT *mt = &r->mt;
239  if (!genrand_initialized(mt)) {
240  r->seed = rand_init(mt, random_seed());
241  }
242  return r;
243 }
244 
245 static struct MT *
247 {
248  return &rand_start(&default_rand)->mt;
249 }
250 
251 unsigned int
253 {
254  struct MT *mt = default_mt();
255  return genrand_int32(mt);
256 }
257 
258 double
260 {
261  struct MT *mt = default_mt();
262  return genrand_real(mt);
263 }
264 
265 #define BDIGITS(x) (RBIGNUM_DIGITS(x))
266 #define BITSPERDIG (SIZEOF_BDIGITS*CHAR_BIT)
267 #define BIGRAD ((BDIGIT_DBL)1 << BITSPERDIG)
268 #define DIGSPERINT (SIZEOF_INT/SIZEOF_BDIGITS)
269 #define BIGUP(x) ((BDIGIT_DBL)(x) << BITSPERDIG)
270 #define BIGDN(x) RSHIFT((x),BITSPERDIG)
271 #define BIGLO(x) ((BDIGIT)((x) & (BIGRAD-1)))
272 #define BDIGMAX ((BDIGIT)-1)
273 
274 #define roomof(n, m) (int)(((n)+(m)-1) / (m))
275 #define numberof(array) (int)(sizeof(array) / sizeof((array)[0]))
276 #define SIZEOF_INT32 (31/CHAR_BIT + 1)
277 
278 static double
279 int_pair_to_real_inclusive(unsigned int a, unsigned int b)
280 {
281  VALUE x = rb_big_new(roomof(64, BITSPERDIG), 1);
282  VALUE m = rb_big_new(roomof(53, BITSPERDIG), 1);
283  BDIGIT *xd = BDIGITS(x);
284  int i = 0;
285  double r;
286 
287  xd[i++] = (BDIGIT)b;
288 #if BITSPERDIG < 32
289  xd[i++] = (BDIGIT)(b >> BITSPERDIG);
290 #endif
291  xd[i++] = (BDIGIT)a;
292 #if BITSPERDIG < 32
293  xd[i++] = (BDIGIT)(a >> BITSPERDIG);
294 #endif
295  xd = BDIGITS(m);
296 #if BITSPERDIG < 53
297  MEMZERO(xd, BDIGIT, roomof(53, BITSPERDIG) - 1);
298 #endif
299  xd[53 / BITSPERDIG] = 1 << 53 % BITSPERDIG;
300  xd[0] |= 1;
301  x = rb_big_mul(x, m);
302  if (FIXNUM_P(x)) {
303 #if CHAR_BIT * SIZEOF_LONG > 64
304  r = (double)(FIX2ULONG(x) >> 64);
305 #else
306  return 0.0;
307 #endif
308  }
309  else {
310 #if 64 % BITSPERDIG == 0
311  long len = RBIGNUM_LEN(x);
312  xd = BDIGITS(x);
313  MEMMOVE(xd, xd + 64 / BITSPERDIG, BDIGIT, len - 64 / BITSPERDIG);
314  MEMZERO(xd + len - 64 / BITSPERDIG, BDIGIT, 64 / BITSPERDIG);
315  r = rb_big2dbl(x);
316 #else
317  x = rb_big_rshift(x, INT2FIX(64));
318  if (FIXNUM_P(x)) {
319  r = (double)FIX2ULONG(x);
320  }
321  else {
322  r = rb_big2dbl(x);
323  }
324 #endif
325  }
326  return ldexp(r, -53);
327 }
328 
331 #define id_minus '-'
332 #define id_plus '+'
334 
335 /* :nodoc: */
336 static void
337 random_mark(void *ptr)
338 {
339  rb_gc_mark(((rb_random_t *)ptr)->seed);
340 }
341 
342 static void
343 random_free(void *ptr)
344 {
345  if (ptr != &default_rand)
346  xfree(ptr);
347 }
348 
349 static size_t
350 random_memsize(const void *ptr)
351 {
352  return ptr ? sizeof(rb_random_t) : 0;
353 }
354 
356  "random",
357  {
358  random_mark,
359  random_free,
361  },
362 };
363 
364 static rb_random_t *
366 {
367  rb_random_t *ptr;
368  TypedData_Get_Struct(obj, rb_random_t, &random_data_type, ptr);
369  return ptr;
370 }
371 
372 static rb_random_t *
374 {
375  if (obj == rb_cRandom) {
376  return rand_start(&default_rand);
377  }
378  if (!rb_typeddata_is_kind_of(obj, &random_data_type)) return NULL;
379  return DATA_PTR(obj);
380 }
381 
382 /* :nodoc: */
383 static VALUE
385 {
386  rb_random_t *rnd;
387  VALUE obj = TypedData_Make_Struct(klass, rb_random_t, &random_data_type, rnd);
388  rnd->seed = INT2FIX(0);
389  return obj;
390 }
391 
392 static VALUE
393 rand_init(struct MT *mt, VALUE vseed)
394 {
395  volatile VALUE seed;
396  long blen = 0;
397  long fixnum_seed;
398  int i, j, len;
399  unsigned int buf0[SIZEOF_LONG / SIZEOF_INT32 * 4], *buf = buf0;
400 
401  seed = rb_to_int(vseed);
402  switch (TYPE(seed)) {
403  case T_FIXNUM:
404  len = 1;
405  fixnum_seed = FIX2LONG(seed);
406  if (fixnum_seed < 0)
407  fixnum_seed = -fixnum_seed;
408  buf[0] = (unsigned int)(fixnum_seed & 0xffffffff);
409 #if SIZEOF_LONG > SIZEOF_INT32
410  if ((long)(int32_t)fixnum_seed != fixnum_seed) {
411  if ((buf[1] = (unsigned int)(fixnum_seed >> 32)) != 0) ++len;
412  }
413 #endif
414  break;
415  case T_BIGNUM:
416  blen = RBIGNUM_LEN(seed);
417  if (blen == 0) {
418  len = 1;
419  }
420  else {
423  len = roomof((int)blen * SIZEOF_BDIGITS, SIZEOF_INT32);
424  }
425  /* allocate ints for init_by_array */
426  if (len > numberof(buf0)) buf = ALLOC_N(unsigned int, len);
427  memset(buf, 0, len * sizeof(*buf));
428  len = 0;
429  for (i = (int)(blen-1); 0 <= i; i--) {
430  j = i * SIZEOF_BDIGITS / SIZEOF_INT32;
431 #if SIZEOF_BDIGITS < SIZEOF_INT32
432  buf[j] <<= BITSPERDIG;
433 #endif
434  buf[j] |= RBIGNUM_DIGITS(seed)[i];
435  if (!len && buf[j]) len = j;
436  }
437  ++len;
438  break;
439  default:
440  rb_raise(rb_eTypeError, "failed to convert %s into Integer",
441  rb_obj_classname(vseed));
442  }
443  if (len <= 1) {
444  init_genrand(mt, buf[0]);
445  }
446  else {
447  if (buf[len-1] == 1) /* remove leading-zero-guard */
448  len--;
449  init_by_array(mt, buf, len);
450  }
451  if (buf != buf0) xfree(buf);
452  return seed;
453 }
454 
455 /*
456  * call-seq: Random.new([seed]) -> prng
457  *
458  * Creates new Mersenne Twister based pseudorandom number generator with
459  * seed. When the argument seed is omitted, the generator is initialized
460  * with Random.new_seed.
461  *
462  * The argument seed is used to ensure repeatable sequences of random numbers
463  * between different runs of the program.
464  *
465  * prng = Random.new(1234)
466  * [ prng.rand, prng.rand ] #=> [0.191519450378892, 0.622108771039832]
467  * [ prng.integer(10), prng.integer(1000) ] #=> [4, 664]
468  * prng = Random.new(1234)
469  * [ prng.rand, prng.rand ] #=> [0.191519450378892, 0.622108771039832]
470  */
471 static VALUE
473 {
474  VALUE vseed;
475  rb_random_t *rnd = get_rnd(obj);
476 
477  if (argc == 0) {
478  vseed = random_seed();
479  }
480  else {
481  rb_scan_args(argc, argv, "01", &vseed);
482  }
483  rnd->seed = rand_init(&rnd->mt, vseed);
484  return obj;
485 }
486 
487 #define DEFAULT_SEED_LEN (DEFAULT_SEED_CNT * (int)sizeof(int))
488 
489 #if defined(S_ISCHR) && !defined(DOSISH)
490 # define USE_DEV_URANDOM 1
491 #else
492 # define USE_DEV_URANDOM 0
493 #endif
494 
495 static void
497 {
498  static int n = 0;
499  struct timeval tv;
500 #if USE_DEV_URANDOM
501  int fd;
502  struct stat statbuf;
503 #elif defined(_WIN32)
504  HCRYPTPROV prov;
505 #endif
506 
507  memset(seed, 0, DEFAULT_SEED_LEN);
508 
509 #if USE_DEV_URANDOM
510  if ((fd = open("/dev/urandom", O_RDONLY
511 #ifdef O_NONBLOCK
512  |O_NONBLOCK
513 #endif
514 #ifdef O_NOCTTY
515  |O_NOCTTY
516 #endif
517  )) >= 0) {
518  rb_update_max_fd(fd);
519  if (fstat(fd, &statbuf) == 0 && S_ISCHR(statbuf.st_mode)) {
520  if (read(fd, seed, DEFAULT_SEED_LEN) < DEFAULT_SEED_LEN) {
521  /* abandon */;
522  }
523  }
524  close(fd);
525  }
526 #elif defined(_WIN32)
527  if (CryptAcquireContext(&prov, NULL, NULL, PROV_RSA_FULL, CRYPT_VERIFYCONTEXT)) {
528  CryptGenRandom(prov, DEFAULT_SEED_LEN, (void *)seed);
529  CryptReleaseContext(prov, 0);
530  }
531 #endif
532 
533  gettimeofday(&tv, 0);
534  seed[0] ^= tv.tv_usec;
535  seed[1] ^= (unsigned int)tv.tv_sec;
536 #if SIZEOF_TIME_T > SIZEOF_INT
537  seed[0] ^= (unsigned int)((time_t)tv.tv_sec >> SIZEOF_INT * CHAR_BIT);
538 #endif
539  seed[2] ^= getpid() ^ (n++ << 16);
540  seed[3] ^= (unsigned int)(VALUE)&seed;
541 #if SIZEOF_VOIDP > SIZEOF_INT
542  seed[2] ^= (unsigned int)((VALUE)&seed >> SIZEOF_INT * CHAR_BIT);
543 #endif
544 }
545 
546 static VALUE
547 make_seed_value(const void *ptr)
548 {
549  const long len = DEFAULT_SEED_LEN/SIZEOF_BDIGITS;
550  BDIGIT *digits;
551  NEWOBJ(big, struct RBignum);
553 
554  RBIGNUM_SET_SIGN(big, 1);
555  rb_big_resize((VALUE)big, len + 1);
556  digits = RBIGNUM_DIGITS(big);
557 
558  MEMCPY(digits, ptr, char, DEFAULT_SEED_LEN);
559 
560  /* set leading-zero-guard if need. */
561  digits[len] =
562 #if SIZEOF_INT32 / SIZEOF_BDIGITS > 1
563  digits[len-2] <= 1 && digits[len-1] == 0
564 #else
565  digits[len-1] <= 1
566 #endif
567  ? 1 : 0;
568 
569  return rb_big_norm((VALUE)big);
570 }
571 
572 /*
573  * call-seq: Random.new_seed -> integer
574  *
575  * Returns arbitrary value for seed.
576  */
577 static VALUE
579 {
580  unsigned int buf[DEFAULT_SEED_CNT];
581  fill_random_seed(buf);
582  return make_seed_value(buf);
583 }
584 
585 /*
586  * call-seq: prng.seed -> integer
587  *
588  * Returns the seed of the generator.
589  */
590 static VALUE
592 {
593  return get_rnd(obj)->seed;
594 }
595 
596 /* :nodoc: */
597 static VALUE
599 {
600  rb_random_t *rnd1 = get_rnd(obj);
601  rb_random_t *rnd2 = get_rnd(orig);
602  struct MT *mt = &rnd1->mt;
603 
604  *rnd1 = *rnd2;
605  mt->next = mt->state + numberof(mt->state) - mt->left + 1;
606  return obj;
607 }
608 
609 static VALUE
610 mt_state(const struct MT *mt)
611 {
612  VALUE bigo = rb_big_new(sizeof(mt->state) / sizeof(BDIGIT), 1);
613  BDIGIT *d = RBIGNUM_DIGITS(bigo);
614  int i;
615 
616  for (i = 0; i < numberof(mt->state); ++i) {
617  unsigned int x = mt->state[i];
618 #if SIZEOF_BDIGITS < SIZEOF_INT32
619  int j;
620  for (j = 0; j < SIZEOF_INT32 / SIZEOF_BDIGITS; ++j) {
621  *d++ = BIGLO(x);
622  x = BIGDN(x);
623  }
624 #else
625  *d++ = (BDIGIT)x;
626 #endif
627  }
628  return rb_big_norm(bigo);
629 }
630 
631 /* :nodoc: */
632 static VALUE
634 {
635  rb_random_t *rnd = get_rnd(obj);
636  return mt_state(&rnd->mt);
637 }
638 
639 /* :nodoc: */
640 static VALUE
642 {
643  return mt_state(&default_rand.mt);
644 }
645 
646 /* :nodoc: */
647 static VALUE
649 {
650  rb_random_t *rnd = get_rnd(obj);
651  return INT2FIX(rnd->mt.left);
652 }
653 
654 /* :nodoc: */
655 static VALUE
657 {
658  return INT2FIX(default_rand.mt.left);
659 }
660 
661 /* :nodoc: */
662 static VALUE
664 {
665  rb_random_t *rnd = get_rnd(obj);
666  VALUE dump = rb_ary_new2(3);
667 
668  rb_ary_push(dump, mt_state(&rnd->mt));
669  rb_ary_push(dump, INT2FIX(rnd->mt.left));
670  rb_ary_push(dump, rnd->seed);
671 
672  return dump;
673 }
674 
675 /* :nodoc: */
676 static VALUE
678 {
679  rb_random_t *rnd = get_rnd(obj);
680  struct MT *mt = &rnd->mt;
681  VALUE state, left = INT2FIX(1), seed = INT2FIX(0);
682  VALUE *ary;
683  unsigned long x;
684 
685  Check_Type(dump, T_ARRAY);
686  ary = RARRAY_PTR(dump);
687  switch (RARRAY_LEN(dump)) {
688  case 3:
689  seed = ary[2];
690  case 2:
691  left = ary[1];
692  case 1:
693  state = ary[0];
694  break;
695  default:
696  rb_raise(rb_eArgError, "wrong dump data");
697  }
698  memset(mt->state, 0, sizeof(mt->state));
699  if (FIXNUM_P(state)) {
700  x = FIX2ULONG(state);
701  mt->state[0] = (unsigned int)x;
702 #if SIZEOF_LONG / SIZEOF_INT >= 2
703  mt->state[1] = (unsigned int)(x >> BITSPERDIG);
704 #endif
705 #if SIZEOF_LONG / SIZEOF_INT >= 3
706  mt->state[2] = (unsigned int)(x >> 2 * BITSPERDIG);
707 #endif
708 #if SIZEOF_LONG / SIZEOF_INT >= 4
709  mt->state[3] = (unsigned int)(x >> 3 * BITSPERDIG);
710 #endif
711  }
712  else {
713  BDIGIT *d;
714  long len;
715  Check_Type(state, T_BIGNUM);
716  len = RBIGNUM_LEN(state);
717  if (len > roomof(sizeof(mt->state), SIZEOF_BDIGITS)) {
718  len = roomof(sizeof(mt->state), SIZEOF_BDIGITS);
719  }
720 #if SIZEOF_BDIGITS < SIZEOF_INT
721  else if (len % DIGSPERINT) {
722  d = RBIGNUM_DIGITS(state) + len;
723 # if DIGSPERINT == 2
724  --len;
725  x = *--d;
726 # else
727  x = 0;
728  do {
729  x = (x << BITSPERDIG) | *--d;
730  } while (--len % DIGSPERINT);
731 # endif
732  mt->state[len / DIGSPERINT] = (unsigned int)x;
733  }
734 #endif
735  if (len > 0) {
736  d = BDIGITS(state) + len;
737  do {
738  --len;
739  x = *--d;
740 # if DIGSPERINT == 2
741  --len;
742  x = (x << BITSPERDIG) | *--d;
743 # elif SIZEOF_BDIGITS < SIZEOF_INT
744  do {
745  x = (x << BITSPERDIG) | *--d;
746  } while (--len % DIGSPERINT);
747 # endif
748  mt->state[len / DIGSPERINT] = (unsigned int)x;
749  } while (len > 0);
750  }
751  }
752  x = NUM2ULONG(left);
753  if (x > numberof(mt->state)) {
754  rb_raise(rb_eArgError, "wrong value");
755  }
756  mt->left = (unsigned int)x;
757  mt->next = mt->state + numberof(mt->state) - x + 1;
758  rnd->seed = rb_to_int(seed);
759 
760  return obj;
761 }
762 
763 /*
764  * call-seq:
765  * srand(number=0) -> old_seed
766  *
767  * Seeds the pseudorandom number generator to the value of
768  * <i>number</i>. If <i>number</i> is omitted,
769  * seeds the generator using a combination of the time, the
770  * process id, and a sequence number. (This is also the behavior if
771  * <code>Kernel::rand</code> is called without previously calling
772  * <code>srand</code>, but without the sequence.) By setting the seed
773  * to a known value, scripts can be made deterministic during testing.
774  * The previous seed value is returned. Also see <code>Kernel::rand</code>.
775  */
776 
777 static VALUE
779 {
780  VALUE seed, old;
782 
783  rb_secure(4);
784  if (argc == 0) {
785  seed = random_seed();
786  }
787  else {
788  rb_scan_args(argc, argv, "01", &seed);
789  }
790  old = r->seed;
791  r->seed = rand_init(&r->mt, seed);
792 
793  return old;
794 }
795 
796 static unsigned long
797 make_mask(unsigned long x)
798 {
799  x = x | x >> 1;
800  x = x | x >> 2;
801  x = x | x >> 4;
802  x = x | x >> 8;
803  x = x | x >> 16;
804 #if 4 < SIZEOF_LONG
805  x = x | x >> 32;
806 #endif
807  return x;
808 }
809 
810 static unsigned long
811 limited_rand(struct MT *mt, unsigned long limit)
812 {
813  /* mt must be initialized */
814  int i;
815  unsigned long val, mask;
816 
817  if (!limit) return 0;
818  mask = make_mask(limit);
819  retry:
820  val = 0;
821  for (i = SIZEOF_LONG/SIZEOF_INT32-1; 0 <= i; i--) {
822  if ((mask >> (i * 32)) & 0xffffffff) {
823  val |= (unsigned long)genrand_int32(mt) << (i * 32);
824  val &= mask;
825  if (limit < val)
826  goto retry;
827  }
828  }
829  return val;
830 }
831 
832 static VALUE
833 limited_big_rand(struct MT *mt, struct RBignum *limit)
834 {
835  /* mt must be initialized */
836  unsigned long mask, lim, rnd;
837  struct RBignum *val;
838  long i, len;
839  int boundary;
840 
841  len = (RBIGNUM_LEN(limit) * SIZEOF_BDIGITS + 3) / 4;
842  val = (struct RBignum *)rb_big_clone((VALUE)limit);
843  RBIGNUM_SET_SIGN(val, 1);
844 #if SIZEOF_BDIGITS == 2
845 # define BIG_GET32(big,i) \
846  (RBIGNUM_DIGITS(big)[(i)*2] | \
847  ((i)*2+1 < RBIGNUM_LEN(big) ? \
848  (RBIGNUM_DIGITS(big)[(i)*2+1] << 16) : \
849  0))
850 # define BIG_SET32(big,i,d) \
851  ((RBIGNUM_DIGITS(big)[(i)*2] = (d) & 0xffff), \
852  ((i)*2+1 < RBIGNUM_LEN(big) ? \
853  (RBIGNUM_DIGITS(big)[(i)*2+1] = (d) >> 16) : \
854  0))
855 #else
856  /* SIZEOF_BDIGITS == 4 */
857 # define BIG_GET32(big,i) (RBIGNUM_DIGITS(big)[(i)])
858 # define BIG_SET32(big,i,d) (RBIGNUM_DIGITS(big)[(i)] = (d))
859 #endif
860  retry:
861  mask = 0;
862  boundary = 1;
863  for (i = len-1; 0 <= i; i--) {
864  lim = BIG_GET32(limit, i);
865  mask = mask ? 0xffffffff : make_mask(lim);
866  if (mask) {
867  rnd = genrand_int32(mt) & mask;
868  if (boundary) {
869  if (lim < rnd)
870  goto retry;
871  if (rnd < lim)
872  boundary = 0;
873  }
874  }
875  else {
876  rnd = 0;
877  }
878  BIG_SET32(val, i, (BDIGIT)rnd);
879  }
880  return rb_big_norm((VALUE)val);
881 }
882 
883 /*
884  * Returns random unsigned long value in [0, _limit_].
885  *
886  * Note that _limit_ is included, and the range of the argument and the
887  * return value depends on environments.
888  */
889 unsigned long
890 rb_genrand_ulong_limited(unsigned long limit)
891 {
892  return limited_rand(default_mt(), limit);
893 }
894 
895 unsigned int
897 {
898  rb_random_t *rnd = try_get_rnd(obj);
899  if (!rnd) {
900 #if SIZEOF_LONG * CHAR_BIT > 32
901  VALUE lim = ULONG2NUM(0x100000000);
902 #elif defined HAVE_LONG_LONG
903  VALUE lim = ULL2NUM((LONG_LONG)0xffffffff+1);
904 #else
905  VALUE lim = rb_big_plus(ULONG2NUM(0xffffffff), INT2FIX(1));
906 #endif
907  return (unsigned int)NUM2ULONG(rb_funcall2(obj, id_rand, 1, &lim));
908  }
909  return genrand_int32(&rnd->mt);
910 }
911 
912 double
914 {
915  rb_random_t *rnd = try_get_rnd(obj);
916  if (!rnd) {
917  VALUE v = rb_funcall2(obj, id_rand, 0, 0);
918  double d = NUM2DBL(v);
919  if (d < 0.0 || d >= 1.0) {
920  rb_raise(rb_eRangeError, "random number too big %g", d);
921  }
922  return d;
923  }
924  return genrand_real(&rnd->mt);
925 }
926 
927 /*
928  * call-seq: prng.bytes(size) -> a_string
929  *
930  * Returns a random binary string. The argument size specified the length of
931  * the result string.
932  */
933 static VALUE
935 {
936  return rb_random_bytes(obj, NUM2LONG(rb_to_int(len)));
937 }
938 
939 VALUE
940 rb_random_bytes(VALUE obj, long n)
941 {
942  rb_random_t *rnd = try_get_rnd(obj);
943  VALUE bytes;
944  char *ptr;
945  unsigned int r, i;
946 
947  if (!rnd) {
948  VALUE len = LONG2NUM(n);
949  return rb_funcall2(obj, id_bytes, 1, &len);
950  }
951  bytes = rb_str_new(0, n);
952  ptr = RSTRING_PTR(bytes);
953  for (; n >= SIZEOF_INT32; n -= SIZEOF_INT32) {
954  r = genrand_int32(&rnd->mt);
955  i = SIZEOF_INT32;
956  do {
957  *ptr++ = (char)r;
958  r >>= CHAR_BIT;
959  } while (--i);
960  }
961  if (n > 0) {
962  r = genrand_int32(&rnd->mt);
963  do {
964  *ptr++ = (char)r;
965  r >>= CHAR_BIT;
966  } while (--n);
967  }
968  return bytes;
969 }
970 
971 static VALUE
972 range_values(VALUE vmax, VALUE *begp, VALUE *endp, int *exclp)
973 {
974  VALUE end, r;
975 
976  if (!rb_range_values(vmax, begp, &end, exclp)) return Qfalse;
977  if (endp) *endp = end;
978  if (!rb_respond_to(end, id_minus)) return Qfalse;
979  r = rb_funcall2(end, id_minus, 1, begp);
980  if (NIL_P(r)) return Qfalse;
981  return r;
982 }
983 
984 static VALUE
985 rand_int(struct MT *mt, VALUE vmax, int restrictive)
986 {
987  /* mt must be initialized */
988  long max;
989  unsigned long r;
990 
991  if (FIXNUM_P(vmax)) {
992  max = FIX2LONG(vmax);
993  if (!max) return Qnil;
994  if (max < 0) {
995  if (restrictive) return Qnil;
996  max = -max;
997  }
998  r = limited_rand(mt, (unsigned long)max - 1);
999  return ULONG2NUM(r);
1000  }
1001  else {
1002  VALUE ret;
1003  if (rb_bigzero_p(vmax)) return Qnil;
1004  if (!RBIGNUM_SIGN(vmax)) {
1005  if (restrictive) return Qnil;
1006  vmax = rb_big_clone(vmax);
1007  RBIGNUM_SET_SIGN(vmax, 1);
1008  }
1009  vmax = rb_big_minus(vmax, INT2FIX(1));
1010  if (FIXNUM_P(vmax)) {
1011  max = FIX2LONG(vmax);
1012  if (max == -1) return Qnil;
1013  r = limited_rand(mt, max);
1014  return LONG2NUM(r);
1015  }
1016  ret = limited_big_rand(mt, RBIGNUM(vmax));
1017  RB_GC_GUARD(vmax);
1018  return ret;
1019  }
1020 }
1021 
1022 static inline double
1024 {
1025  double x = RFLOAT_VALUE(v);
1026  if (isinf(x) || isnan(x)) {
1027  VALUE error = INT2FIX(EDOM);
1029  }
1030  return x;
1031 }
1032 
1033 static inline VALUE
1034 rand_range(struct MT* mt, VALUE range)
1035 {
1036  VALUE beg = Qundef, end = Qundef, vmax, v;
1037  int excl = 0;
1038 
1039  if ((v = vmax = range_values(range, &beg, &end, &excl)) == Qfalse)
1040  return Qfalse;
1041  if (TYPE(vmax) != T_FLOAT && (v = rb_check_to_integer(vmax, "to_int"), !NIL_P(v))) {
1042  long max;
1043  vmax = v;
1044  v = Qnil;
1045  if (FIXNUM_P(vmax)) {
1046  fixnum:
1047  if ((max = FIX2LONG(vmax) - excl) >= 0) {
1048  unsigned long r = limited_rand(mt, (unsigned long)max);
1049  v = ULONG2NUM(r);
1050  }
1051  }
1052  else if (BUILTIN_TYPE(vmax) == T_BIGNUM && RBIGNUM_SIGN(vmax) && !rb_bigzero_p(vmax)) {
1053  vmax = excl ? rb_big_minus(vmax, INT2FIX(1)) : rb_big_norm(vmax);
1054  if (FIXNUM_P(vmax)) {
1055  excl = 0;
1056  goto fixnum;
1057  }
1058  v = limited_big_rand(mt, RBIGNUM(vmax));
1059  }
1060  }
1061  else if (v = rb_check_to_float(vmax), !NIL_P(v)) {
1062  int scale = 1;
1063  double max = RFLOAT_VALUE(v), mid = 0.5, r;
1064  if (isinf(max)) {
1065  double min = float_value(rb_to_float(beg)) / 2.0;
1066  max = float_value(rb_to_float(end)) / 2.0;
1067  scale = 2;
1068  mid = max + min;
1069  max -= min;
1070  }
1071  else {
1072  float_value(v);
1073  }
1074  v = Qnil;
1075  if (max > 0.0) {
1076  if (excl) {
1077  r = genrand_real(mt);
1078  }
1079  else {
1080  r = genrand_real2(mt);
1081  }
1082  if (scale > 1) {
1083  return rb_float_new(+(+(+(r - 0.5) * max) * scale) + mid);
1084  }
1085  v = rb_float_new(r * max);
1086  }
1087  else if (max == 0.0 && !excl) {
1088  v = rb_float_new(0.0);
1089  }
1090  }
1091 
1092  if (FIXNUM_P(beg) && FIXNUM_P(v)) {
1093  long x = FIX2LONG(beg) + FIX2LONG(v);
1094  return LONG2NUM(x);
1095  }
1096  switch (TYPE(v)) {
1097  case T_NIL:
1098  break;
1099  case T_BIGNUM:
1100  return rb_big_plus(v, beg);
1101  case T_FLOAT: {
1102  VALUE f = rb_check_to_float(beg);
1103  if (!NIL_P(f)) {
1104  RFLOAT_VALUE(v) += RFLOAT_VALUE(f);
1105  return v;
1106  }
1107  }
1108  default:
1109  return rb_funcall2(beg, id_plus, 1, &v);
1110  }
1111 
1112  return v;
1113 }
1114 
1115 /*
1116  * call-seq:
1117  * prng.rand -> float
1118  * prng.rand(limit) -> number
1119  *
1120  * When the argument is an +Integer+ or a +Bignum+, it returns a
1121  * random integer greater than or equal to zero and less than the
1122  * argument. Unlike Random.rand, when the argument is a negative
1123  * integer or zero, it raises an ArgumentError.
1124  *
1125  * When the argument is a +Float+, it returns a random floating point
1126  * number between 0.0 and _max_, including 0.0 and excluding _max_.
1127  *
1128  * When the argument _limit_ is a +Range+, it returns a random
1129  * number where range.member?(number) == true.
1130  * prng.rand(5..9) #=> one of [5, 6, 7, 8, 9]
1131  * prng.rand(5...9) #=> one of [5, 6, 7, 8]
1132  * prng.rand(5.0..9.0) #=> between 5.0 and 9.0, including 9.0
1133  * prng.rand(5.0...9.0) #=> between 5.0 and 9.0, excluding 9.0
1134  *
1135  * +begin+/+end+ of the range have to have subtract and add methods.
1136  *
1137  * Otherwise, it raises an ArgumentError.
1138  */
1139 static VALUE
1141 {
1142  rb_random_t *rnd = get_rnd(obj);
1143  VALUE vmax, v;
1144 
1145  if (argc == 0) {
1146  return rb_float_new(genrand_real(&rnd->mt));
1147  }
1148  else if (argc != 1) {
1149  rb_raise(rb_eArgError, "wrong number of arguments (%d for 0..1)", argc);
1150  }
1151  vmax = argv[0];
1152  if (NIL_P(vmax)) {
1153  v = Qnil;
1154  }
1155  else if (TYPE(vmax) != T_FLOAT && (v = rb_check_to_integer(vmax, "to_int"), !NIL_P(v))) {
1156  v = rand_int(&rnd->mt, v, 1);
1157  }
1158  else if (v = rb_check_to_float(vmax), !NIL_P(v)) {
1159  double max = float_value(v);
1160  if (max > 0.0)
1161  v = rb_float_new(max * genrand_real(&rnd->mt));
1162  else
1163  v = Qnil;
1164  }
1165  else if ((v = rand_range(&rnd->mt, vmax)) != Qfalse) {
1166  /* nothing to do */
1167  }
1168  else {
1169  v = Qnil;
1170  (void)NUM2LONG(vmax);
1171  }
1172  if (NIL_P(v)) {
1173  VALUE mesg = rb_str_new_cstr("invalid argument - ");
1174  rb_str_append(mesg, rb_obj_as_string(argv[0]));
1176  }
1177 
1178  return v;
1179 }
1180 
1181 /*
1182  * call-seq:
1183  * prng1 == prng2 -> true or false
1184  *
1185  * Returns true if the generators' states equal.
1186  */
1187 static VALUE
1189 {
1190  rb_random_t *r1, *r2;
1191  if (rb_obj_class(self) != rb_obj_class(other)) return Qfalse;
1192  r1 = get_rnd(self);
1193  r2 = get_rnd(other);
1194  if (!RTEST(rb_funcall2(r1->seed, rb_intern("=="), 1, &r2->seed))) return Qfalse;
1195  if (memcmp(r1->mt.state, r2->mt.state, sizeof(r1->mt.state))) return Qfalse;
1196  if ((r1->mt.next - r1->mt.state) != (r2->mt.next - r2->mt.state)) return Qfalse;
1197  if (r1->mt.left != r2->mt.left) return Qfalse;
1198  return Qtrue;
1199 }
1200 
1201 /*
1202  * call-seq:
1203  * rand(max=0) -> number
1204  *
1205  *
1206  * If <i>max</i> is +Range+, returns a pseudorandom number where
1207  * range.member(number) == true.
1208  *
1209  * Or else converts _max_ to an integer using max1 =
1210  * max<code>.to_i.abs</code>.
1211  *
1212  * Then if _max_ is +nil+ the result is zero, returns a pseudorandom floating
1213  * point number greater than or equal to 0.0 and less than 1.0.
1214  *
1215  * Otherwise, returns a pseudorandom integer greater than or equal to zero and
1216  * less than max1.
1217  *
1218  * <code>Kernel::srand</code> may be used to ensure repeatable sequences of
1219  * random numbers between different runs of the program. Ruby currently uses
1220  * a modified Mersenne Twister with a period of 2**19937-1.
1221  *
1222  * srand 1234 #=> 0
1223  * [ rand, rand ] #=> [0.191519450163469, 0.49766366626136]
1224  * [ rand(10), rand(1000) ] #=> [6, 817]
1225  * srand 1234 #=> 1234
1226  * [ rand, rand ] #=> [0.191519450163469, 0.49766366626136]
1227  */
1228 
1229 static VALUE
1231 {
1232  VALUE v, vmax, r;
1233  struct MT *mt = default_mt();
1234 
1235  if (argc == 0) goto zero_arg;
1236  rb_scan_args(argc, argv, "01", &vmax);
1237  if (NIL_P(vmax)) goto zero_arg;
1238  if ((v = rand_range(mt, vmax)) != Qfalse) {
1239  return v;
1240  }
1241  vmax = rb_to_int(vmax);
1242  if (vmax == INT2FIX(0) || NIL_P(r = rand_int(mt, vmax, 0))) {
1243  zero_arg:
1244  return DBL2NUM(genrand_real(mt));
1245  }
1246  return r;
1247 }
1248 
1249 /*
1250  * call-seq:
1251  * Random.rand -> float
1252  * Random.rand(limit) -> number
1253  *
1254  * Alias of _Random::DEFAULT.rand_.
1255  *
1256  */
1257 
1258 static VALUE
1260 {
1261  rand_start(&default_rand);
1262  return random_rand(argc, argv, rb_Random_DEFAULT);
1263 }
1264 
1265 #define SIP_HASH_STREAMING 0
1266 #define sip_hash24 ruby_sip_hash24
1267 #if !defined _WIN32 && !defined BYTE_ORDER
1268 # ifdef WORDS_BIGENDIAN
1269 # define BYTE_ORDER BIG_ENDIAN
1270 # else
1271 # define BYTE_ORDER LITTLE_ENDIAN
1272 # endif
1273 # ifndef LITTLE_ENDIAN
1274 # define LITTLE_ENDIAN 1234
1275 # endif
1276 # ifndef BIG_ENDIAN
1277 # define BIG_ENDIAN 4321
1278 # endif
1279 #endif
1280 #include "siphash.c"
1281 
1283 static union {
1285  uint32_t u32[(16 * sizeof(uint8_t) - 1) / sizeof(uint32_t)];
1286 } sipseed;
1287 
1288 static VALUE
1289 init_randomseed(struct MT *mt, unsigned int initial[DEFAULT_SEED_CNT])
1290 {
1291  VALUE seed;
1292  fill_random_seed(initial);
1293  init_by_array(mt, initial, DEFAULT_SEED_CNT);
1294  seed = make_seed_value(initial);
1295  memset(initial, 0, DEFAULT_SEED_LEN);
1296  return seed;
1297 }
1298 
1299 void
1301 {
1302  rb_random_t *r = &default_rand;
1303  unsigned int initial[DEFAULT_SEED_CNT];
1304  struct MT *mt = &r->mt;
1305  VALUE seed = init_randomseed(mt, initial);
1306  int i;
1307 
1308  hashseed = genrand_int32(mt);
1309 #if SIZEOF_ST_INDEX_T*CHAR_BIT > 4*8
1310  hashseed <<= 32;
1311  hashseed |= genrand_int32(mt);
1312 #endif
1313 #if SIZEOF_ST_INDEX_T*CHAR_BIT > 8*8
1314  hashseed <<= 32;
1315  hashseed |= genrand_int32(mt);
1316 #endif
1317 #if SIZEOF_ST_INDEX_T*CHAR_BIT > 12*8
1318  hashseed <<= 32;
1319  hashseed |= genrand_int32(mt);
1320 #endif
1321 
1322  for (i = 0; i < numberof(sipseed.u32); ++i)
1323  sipseed.u32[i] = genrand_int32(mt);
1324 
1325  rb_global_variable(&r->seed);
1326  r->seed = seed;
1327 }
1328 
1329 st_index_t
1331 {
1332  return st_hash_start(hashseed + h);
1333 }
1334 
1335 st_index_t
1336 rb_memhash(const void *ptr, long len)
1337 {
1338  sip_uint64_t h = sip_hash24(sipseed.key, ptr, len);
1339 #ifdef HAVE_UINT64_T
1340  return (st_index_t)h;
1341 #else
1342  return (st_index_t)(h.u32[0] ^ h.u32[1]);
1343 #endif
1344 }
1345 
1346 static void
1348 {
1349  VALUE seed = default_rand.seed;
1350 
1351  if (RB_TYPE_P(seed, T_BIGNUM)) {
1352  RBASIC(seed)->klass = rb_cBignum;
1353  }
1354 }
1355 
1356 void
1358 {
1359  rb_random_t *r = &default_rand;
1360  uninit_genrand(&r->mt);
1361  r->seed = INT2FIX(0);
1362 }
1363 
1364 void
1366 {
1367  Init_RandomSeed2();
1368  rb_define_global_function("srand", rb_f_srand, -1);
1369  rb_define_global_function("rand", rb_f_rand, -1);
1370 
1371  rb_cRandom = rb_define_class("Random", rb_cObject);
1373  rb_define_method(rb_cRandom, "initialize", random_init, -1);
1374  rb_define_method(rb_cRandom, "rand", random_rand, -1);
1377  rb_define_method(rb_cRandom, "initialize_copy", random_copy, 1);
1378  rb_define_method(rb_cRandom, "marshal_dump", random_dump, 0);
1379  rb_define_method(rb_cRandom, "marshal_load", random_load, 1);
1383 
1384  rb_Random_DEFAULT = TypedData_Wrap_Struct(rb_cRandom, &random_data_type, &default_rand);
1387 
1393 
1394  id_rand = rb_intern("rand");
1395  id_bytes = rb_intern("bytes");
1396 }
1397