1 /* Random objects */
2
3 /* ------------------------------------------------------------------
4 The code in this module was based on a download from:
5 http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html
6
7 It was modified in 2002 by Raymond Hettinger as follows:
8
9 * the principal computational lines untouched.
10
11 * renamed genrand_res53() to random_random() and wrapped
12 in python calling/return code.
13
14 * genrand_int32() and the helper functions, init_genrand()
15 and init_by_array(), were declared static, wrapped in
16 Python calling/return code. also, their global data
17 references were replaced with structure references.
18
19 * unused functions from the original were deleted.
20 new, original C python code was added to implement the
21 Random() interface.
22
23 The following are the verbatim comments from the original code:
24
25 A C-program for MT19937, with initialization improved 2002/1/26.
26 Coded by Takuji Nishimura and Makoto Matsumoto.
27
28 Before using, initialize the state by using init_genrand(seed)
29 or init_by_array(init_key, key_length).
30
31 Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
32 All rights reserved.
33
34 Redistribution and use in source and binary forms, with or without
35 modification, are permitted provided that the following conditions
36 are met:
37
38 1. Redistributions of source code must retain the above copyright
39 notice, this list of conditions and the following disclaimer.
40
41 2. Redistributions in binary form must reproduce the above copyright
42 notice, this list of conditions and the following disclaimer in the
43 documentation and/or other materials provided with the distribution.
44
45 3. The names of its contributors may not be used to endorse or promote
46 products derived from this software without specific prior written
47 permission.
48
49 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
50 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
51 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
52 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
53 CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
54 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
55 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
56 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
57 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
58 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
59 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
60
61
62 Any feedback is very welcome.
63 http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
64 email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
65 */
66
67 /* ---------------------------------------------------------------*/
68
69 #include "Python.h"
70 #include <time.h> /* for seeding to current time */
71
72 /* Period parameters -- These are all magic. Don't change. */
73 #define N 624
74 #define M 397
75 #define MATRIX_A 0x9908b0dfUL /* constant vector a */
76 #define UPPER_MASK 0x80000000UL /* most significant w-r bits */
77 #define LOWER_MASK 0x7fffffffUL /* least significant r bits */
78
79 typedef struct {
80 PyObject_HEAD
81 unsigned long state[N];
82 int index;
83 } RandomObject;
84
85 static PyTypeObject Random_Type;
86
87 #define RandomObject_Check(v) (Py_TYPE(v) == &Random_Type)
88
89
90 /* Random methods */
91
92
93 /* generates a random number on [0,0xffffffff]-interval */
94 static unsigned long
genrand_int32(RandomObject * self)95 genrand_int32(RandomObject *self)
96 {
97 unsigned long y;
98 static unsigned long mag01[2]={0x0UL, MATRIX_A};
99 /* mag01[x] = x * MATRIX_A for x=0,1 */
100 unsigned long *mt;
101
102 mt = self->state;
103 if (self->index >= N) { /* generate N words at one time */
104 int kk;
105
106 for (kk=0;kk<N-M;kk++) {
107 y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
108 mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];
109 }
110 for (;kk<N-1;kk++) {
111 y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
112 mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
113 }
114 y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
115 mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
116
117 self->index = 0;
118 }
119
120 y = mt[self->index++];
121 y ^= (y >> 11);
122 y ^= (y << 7) & 0x9d2c5680UL;
123 y ^= (y << 15) & 0xefc60000UL;
124 y ^= (y >> 18);
125 return y;
126 }
127
128 /* random_random is the function named genrand_res53 in the original code;
129 * generates a random number on [0,1) with 53-bit resolution; note that
130 * 9007199254740992 == 2**53; I assume they're spelling "/2**53" as
131 * multiply-by-reciprocal in the (likely vain) hope that the compiler will
132 * optimize the division away at compile-time. 67108864 is 2**26. In
133 * effect, a contains 27 random bits shifted left 26, and b fills in the
134 * lower 26 bits of the 53-bit numerator.
135 * The orginal code credited Isaku Wada for this algorithm, 2002/01/09.
136 */
137 static PyObject *
random_random(RandomObject * self)138 random_random(RandomObject *self)
139 {
140 unsigned long a=genrand_int32(self)>>5, b=genrand_int32(self)>>6;
141 return PyFloat_FromDouble((a*67108864.0+b)*(1.0/9007199254740992.0));
142 }
143
144 /* initializes mt[N] with a seed */
145 static void
init_genrand(RandomObject * self,unsigned long s)146 init_genrand(RandomObject *self, unsigned long s)
147 {
148 int mti;
149 unsigned long *mt;
150
151 mt = self->state;
152 mt[0]= s & 0xffffffffUL;
153 for (mti=1; mti<N; mti++) {
154 mt[mti] =
155 (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
156 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
157 /* In the previous versions, MSBs of the seed affect */
158 /* only MSBs of the array mt[]. */
159 /* 2002/01/09 modified by Makoto Matsumoto */
160 mt[mti] &= 0xffffffffUL;
161 /* for >32 bit machines */
162 }
163 self->index = mti;
164 return;
165 }
166
167 /* initialize by an array with array-length */
168 /* init_key is the array for initializing keys */
169 /* key_length is its length */
170 static PyObject *
init_by_array(RandomObject * self,unsigned long init_key[],unsigned long key_length)171 init_by_array(RandomObject *self, unsigned long init_key[], unsigned long key_length)
172 {
173 unsigned int i, j, k; /* was signed in the original code. RDH 12/16/2002 */
174 unsigned long *mt;
175
176 mt = self->state;
177 init_genrand(self, 19650218UL);
178 i=1; j=0;
179 k = (N>key_length ? N : key_length);
180 for (; k; k--) {
181 mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL))
182 + init_key[j] + j; /* non linear */
183 mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
184 i++; j++;
185 if (i>=N) { mt[0] = mt[N-1]; i=1; }
186 if (j>=key_length) j=0;
187 }
188 for (k=N-1; k; k--) {
189 mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL))
190 - i; /* non linear */
191 mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
192 i++;
193 if (i>=N) { mt[0] = mt[N-1]; i=1; }
194 }
195
196 mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
197 Py_INCREF(Py_None);
198 return Py_None;
199 }
200
201 /*
202 * The rest is Python-specific code, neither part of, nor derived from, the
203 * Twister download.
204 */
205
206 static PyObject *
random_seed(RandomObject * self,PyObject * args)207 random_seed(RandomObject *self, PyObject *args)
208 {
209 PyObject *result = NULL; /* guilty until proved innocent */
210 PyObject *masklower = NULL;
211 PyObject *thirtytwo = NULL;
212 PyObject *n = NULL;
213 unsigned long *key = NULL;
214 unsigned long keymax; /* # of allocated slots in key */
215 unsigned long keyused; /* # of used slots in key */
216 int err;
217
218 PyObject *arg = NULL;
219
220 if (!PyArg_UnpackTuple(args, "seed", 0, 1, &arg))
221 return NULL;
222
223 if (arg == NULL || arg == Py_None) {
224 time_t now;
225
226 time(&now);
227 init_genrand(self, (unsigned long)now);
228 Py_INCREF(Py_None);
229 return Py_None;
230 }
231 /* If the arg is an int or long, use its absolute value; else use
232 * the absolute value of its hash code.
233 * Calling int.__abs__() or long.__abs__() prevents calling arg.__abs__(),
234 * which might return an invalid value. See issue #31478.
235 */
236 if (PyInt_Check(arg)) {
237 n = PyInt_Type.tp_as_number->nb_absolute(arg);
238 }
239 else if (PyLong_Check(arg)) {
240 n = PyLong_Type.tp_as_number->nb_absolute(arg);
241 }
242 else {
243 long hash = PyObject_Hash(arg);
244 if (hash == -1)
245 goto Done;
246 n = PyLong_FromUnsignedLong((unsigned long)hash);
247 }
248 if (n == NULL)
249 goto Done;
250
251 /* Now split n into 32-bit chunks, from the right. Each piece is
252 * stored into key, which has a capacity of keymax chunks, of which
253 * keyused are filled. Alas, the repeated shifting makes this a
254 * quadratic-time algorithm; we'd really like to use
255 * _PyLong_AsByteArray here, but then we'd have to break into the
256 * long representation to figure out how big an array was needed
257 * in advance.
258 */
259 keymax = 8; /* arbitrary; grows later if needed */
260 keyused = 0;
261 key = (unsigned long *)PyMem_Malloc(keymax * sizeof(*key));
262 if (key == NULL)
263 goto Done;
264
265 masklower = PyLong_FromUnsignedLong(0xffffffffU);
266 if (masklower == NULL)
267 goto Done;
268 thirtytwo = PyInt_FromLong(32L);
269 if (thirtytwo == NULL)
270 goto Done;
271 while ((err=PyObject_IsTrue(n))) {
272 PyObject *newn;
273 PyObject *pychunk;
274 unsigned long chunk;
275
276 if (err == -1)
277 goto Done;
278 pychunk = PyNumber_And(n, masklower);
279 if (pychunk == NULL)
280 goto Done;
281 chunk = PyLong_AsUnsignedLong(pychunk);
282 Py_DECREF(pychunk);
283 if (chunk == (unsigned long)-1 && PyErr_Occurred())
284 goto Done;
285 newn = PyNumber_Rshift(n, thirtytwo);
286 if (newn == NULL)
287 goto Done;
288 Py_DECREF(n);
289 n = newn;
290 if (keyused >= keymax) {
291 unsigned long bigger = keymax << 1;
292 if ((bigger >> 1) != keymax) {
293 PyErr_NoMemory();
294 goto Done;
295 }
296 key = (unsigned long *)PyMem_Realloc(key,
297 bigger * sizeof(*key));
298 if (key == NULL)
299 goto Done;
300 keymax = bigger;
301 }
302 assert(keyused < keymax);
303 key[keyused++] = chunk;
304 }
305
306 if (keyused == 0)
307 key[keyused++] = 0UL;
308 result = init_by_array(self, key, keyused);
309 Done:
310 Py_XDECREF(masklower);
311 Py_XDECREF(thirtytwo);
312 Py_XDECREF(n);
313 PyMem_Free(key);
314 return result;
315 }
316
317 static PyObject *
random_getstate(RandomObject * self)318 random_getstate(RandomObject *self)
319 {
320 PyObject *state;
321 PyObject *element;
322 int i;
323
324 state = PyTuple_New(N+1);
325 if (state == NULL)
326 return NULL;
327 for (i=0; i<N ; i++) {
328 element = PyLong_FromUnsignedLong(self->state[i]);
329 if (element == NULL)
330 goto Fail;
331 PyTuple_SET_ITEM(state, i, element);
332 }
333 element = PyLong_FromLong((long)(self->index));
334 if (element == NULL)
335 goto Fail;
336 PyTuple_SET_ITEM(state, i, element);
337 return state;
338
339 Fail:
340 Py_DECREF(state);
341 return NULL;
342 }
343
344 static PyObject *
random_setstate(RandomObject * self,PyObject * state)345 random_setstate(RandomObject *self, PyObject *state)
346 {
347 int i;
348 unsigned long element;
349 long index;
350 unsigned long new_state[N];
351
352 if (!PyTuple_Check(state)) {
353 PyErr_SetString(PyExc_TypeError,
354 "state vector must be a tuple");
355 return NULL;
356 }
357 if (PyTuple_Size(state) != N+1) {
358 PyErr_SetString(PyExc_ValueError,
359 "state vector is the wrong size");
360 return NULL;
361 }
362
363 for (i=0; i<N ; i++) {
364 element = PyLong_AsUnsignedLong(PyTuple_GET_ITEM(state, i));
365 if (element == (unsigned long)-1 && PyErr_Occurred())
366 return NULL;
367 new_state[i] = element & 0xffffffffUL; /* Make sure we get sane state */
368 }
369
370 index = PyLong_AsLong(PyTuple_GET_ITEM(state, i));
371 if (index == -1 && PyErr_Occurred())
372 return NULL;
373 if (index < 0 || index > N) {
374 PyErr_SetString(PyExc_ValueError, "invalid state");
375 return NULL;
376 }
377 self->index = (int)index;
378 for (i = 0; i < N; i++)
379 self->state[i] = new_state[i];
380
381 Py_INCREF(Py_None);
382 return Py_None;
383 }
384
385 /*
386 Jumpahead should be a fast way advance the generator n-steps ahead, but
387 lacking a formula for that, the next best is to use n and the existing
388 state to create a new state far away from the original.
389
390 The generator uses constant spaced additive feedback, so shuffling the
391 state elements ought to produce a state which would not be encountered
392 (in the near term) by calls to random(). Shuffling is normally
393 implemented by swapping the ith element with another element ranging
394 from 0 to i inclusive. That allows the element to have the possibility
395 of not being moved. Since the goal is to produce a new, different
396 state, the swap element is ranged from 0 to i-1 inclusive. This assures
397 that each element gets moved at least once.
398
399 To make sure that consecutive calls to jumpahead(n) produce different
400 states (even in the rare case of involutory shuffles), i+1 is added to
401 each element at position i. Successive calls are then guaranteed to
402 have changing (growing) values as well as shuffled positions.
403
404 Finally, the self->index value is set to N so that the generator itself
405 kicks in on the next call to random(). This assures that all results
406 have been through the generator and do not just reflect alterations to
407 the underlying state.
408 */
409
410 static PyObject *
random_jumpahead(RandomObject * self,PyObject * n)411 random_jumpahead(RandomObject *self, PyObject *n)
412 {
413 long i, j;
414 PyObject *iobj;
415 PyObject *remobj;
416 unsigned long *mt, tmp, nonzero;
417
418 if (!_PyAnyInt_Check(n)) {
419 PyErr_Format(PyExc_TypeError, "jumpahead requires an "
420 "integer, not '%s'",
421 Py_TYPE(n)->tp_name);
422 return NULL;
423 }
424
425 mt = self->state;
426 for (i = N-1; i > 1; i--) {
427 iobj = PyInt_FromLong(i);
428 if (iobj == NULL)
429 return NULL;
430 remobj = PyNumber_Remainder(n, iobj);
431 Py_DECREF(iobj);
432 if (remobj == NULL)
433 return NULL;
434 j = PyInt_AsLong(remobj);
435 Py_DECREF(remobj);
436 if (j == -1L && PyErr_Occurred())
437 return NULL;
438 tmp = mt[i];
439 mt[i] = mt[j];
440 mt[j] = tmp;
441 }
442
443 nonzero = 0;
444 for (i = 1; i < N; i++) {
445 mt[i] += i+1;
446 mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
447 nonzero |= mt[i];
448 }
449
450 /* Ensure the state is nonzero: in the unlikely event that mt[1] through
451 mt[N-1] are all zero, set the MSB of mt[0] (see issue #14591). In the
452 normal case, we fall back to the pre-issue 14591 behaviour for mt[0]. */
453 if (nonzero) {
454 mt[0] += 1;
455 mt[0] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
456 }
457 else {
458 mt[0] = 0x80000000UL;
459 }
460
461 self->index = N;
462 Py_INCREF(Py_None);
463 return Py_None;
464 }
465
466 static PyObject *
random_getrandbits(RandomObject * self,PyObject * args)467 random_getrandbits(RandomObject *self, PyObject *args)
468 {
469 int k, i, bytes;
470 unsigned long r;
471 unsigned char *bytearray;
472 PyObject *result;
473
474 if (!PyArg_ParseTuple(args, "i:getrandbits", &k))
475 return NULL;
476
477 if (k <= 0) {
478 PyErr_SetString(PyExc_ValueError,
479 "number of bits must be greater than zero");
480 return NULL;
481 }
482
483 bytes = ((k - 1) / 32 + 1) * 4;
484 bytearray = (unsigned char *)PyMem_Malloc(bytes);
485 if (bytearray == NULL) {
486 PyErr_NoMemory();
487 return NULL;
488 }
489
490 /* Fill-out whole words, byte-by-byte to avoid endianness issues */
491 for (i=0 ; i<bytes ; i+=4, k-=32) {
492 r = genrand_int32(self);
493 if (k < 32)
494 r >>= (32 - k);
495 bytearray[i+0] = (unsigned char)r;
496 bytearray[i+1] = (unsigned char)(r >> 8);
497 bytearray[i+2] = (unsigned char)(r >> 16);
498 bytearray[i+3] = (unsigned char)(r >> 24);
499 }
500
501 /* little endian order to match bytearray assignment order */
502 result = _PyLong_FromByteArray(bytearray, bytes, 1, 0);
503 PyMem_Free(bytearray);
504 return result;
505 }
506
507 static PyObject *
random_new(PyTypeObject * type,PyObject * args,PyObject * kwds)508 random_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
509 {
510 RandomObject *self;
511 PyObject *tmp;
512
513 if (type == &Random_Type && !_PyArg_NoKeywords("Random()", kwds))
514 return NULL;
515
516 self = (RandomObject *)type->tp_alloc(type, 0);
517 if (self == NULL)
518 return NULL;
519 tmp = random_seed(self, args);
520 if (tmp == NULL) {
521 Py_DECREF(self);
522 return NULL;
523 }
524 Py_DECREF(tmp);
525 return (PyObject *)self;
526 }
527
528 static PyMethodDef random_methods[] = {
529 {"random", (PyCFunction)random_random, METH_NOARGS,
530 PyDoc_STR("random() -> x in the interval [0, 1).")},
531 {"seed", (PyCFunction)random_seed, METH_VARARGS,
532 PyDoc_STR("seed([n]) -> None. Defaults to current time.")},
533 {"getstate", (PyCFunction)random_getstate, METH_NOARGS,
534 PyDoc_STR("getstate() -> tuple containing the current state.")},
535 {"setstate", (PyCFunction)random_setstate, METH_O,
536 PyDoc_STR("setstate(state) -> None. Restores generator state.")},
537 {"jumpahead", (PyCFunction)random_jumpahead, METH_O,
538 PyDoc_STR("jumpahead(int) -> None. Create new state from "
539 "existing state and integer.")},
540 {"getrandbits", (PyCFunction)random_getrandbits, METH_VARARGS,
541 PyDoc_STR("getrandbits(k) -> x. Generates a long int with "
542 "k random bits.")},
543 {NULL, NULL} /* sentinel */
544 };
545
546 PyDoc_STRVAR(random_doc,
547 "Random() -> create a random number generator with its own internal state.");
548
549 static PyTypeObject Random_Type = {
550 PyVarObject_HEAD_INIT(NULL, 0)
551 "_random.Random", /*tp_name*/
552 sizeof(RandomObject), /*tp_basicsize*/
553 0, /*tp_itemsize*/
554 /* methods */
555 0, /*tp_dealloc*/
556 0, /*tp_print*/
557 0, /*tp_getattr*/
558 0, /*tp_setattr*/
559 0, /*tp_compare*/
560 0, /*tp_repr*/
561 0, /*tp_as_number*/
562 0, /*tp_as_sequence*/
563 0, /*tp_as_mapping*/
564 0, /*tp_hash*/
565 0, /*tp_call*/
566 0, /*tp_str*/
567 PyObject_GenericGetAttr, /*tp_getattro*/
568 0, /*tp_setattro*/
569 0, /*tp_as_buffer*/
570 Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /*tp_flags*/
571 random_doc, /*tp_doc*/
572 0, /*tp_traverse*/
573 0, /*tp_clear*/
574 0, /*tp_richcompare*/
575 0, /*tp_weaklistoffset*/
576 0, /*tp_iter*/
577 0, /*tp_iternext*/
578 random_methods, /*tp_methods*/
579 0, /*tp_members*/
580 0, /*tp_getset*/
581 0, /*tp_base*/
582 0, /*tp_dict*/
583 0, /*tp_descr_get*/
584 0, /*tp_descr_set*/
585 0, /*tp_dictoffset*/
586 0, /*tp_init*/
587 0, /*tp_alloc*/
588 random_new, /*tp_new*/
589 _PyObject_Del, /*tp_free*/
590 0, /*tp_is_gc*/
591 };
592
593 PyDoc_STRVAR(module_doc,
594 "Module implements the Mersenne Twister random number generator.");
595
596 PyMODINIT_FUNC
init_random(void)597 init_random(void)
598 {
599 PyObject *m;
600
601 if (PyType_Ready(&Random_Type) < 0)
602 return;
603 m = Py_InitModule3("_random", NULL, module_doc);
604 if (m == NULL)
605 return;
606 Py_INCREF(&Random_Type);
607 PyModule_AddObject(m, "Random", (PyObject *)&Random_Type);
608 }
609