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