Skip to content

HTTPS clone URL

Subversion checkout URL

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