No issues found
1 /* Random objects */
2
3 /* ------------------------------------------------------------------
4 The code in this module was based on a download from:
5 http://www.math.keio.ac.jp/~matumoto/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.keio.ac.jp/matumoto/emt.html
64 email: matumoto@math.keio.ac.jp
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
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 *
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
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 *
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 *
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 */
234 if (PyInt_Check(arg) || PyLong_Check(arg))
235 n = PyNumber_Absolute(arg);
236 else {
237 long hash = PyObject_Hash(arg);
238 if (hash == -1)
239 goto Done;
240 n = PyLong_FromUnsignedLong((unsigned long)hash);
241 }
242 if (n == NULL)
243 goto Done;
244
245 /* Now split n into 32-bit chunks, from the right. Each piece is
246 * stored into key, which has a capacity of keymax chunks, of which
247 * keyused are filled. Alas, the repeated shifting makes this a
248 * quadratic-time algorithm; we'd really like to use
249 * _PyLong_AsByteArray here, but then we'd have to break into the
250 * long representation to figure out how big an array was needed
251 * in advance.
252 */
253 keymax = 8; /* arbitrary; grows later if needed */
254 keyused = 0;
255 key = (unsigned long *)PyMem_Malloc(keymax * sizeof(*key));
256 if (key == NULL)
257 goto Done;
258
259 masklower = PyLong_FromUnsignedLong(0xffffffffU);
260 if (masklower == NULL)
261 goto Done;
262 thirtytwo = PyInt_FromLong(32L);
263 if (thirtytwo == NULL)
264 goto Done;
265 while ((err=PyObject_IsTrue(n))) {
266 PyObject *newn;
267 PyObject *pychunk;
268 unsigned long chunk;
269
270 if (err == -1)
271 goto Done;
272 pychunk = PyNumber_And(n, masklower);
273 if (pychunk == NULL)
274 goto Done;
275 chunk = PyLong_AsUnsignedLong(pychunk);
276 Py_DECREF(pychunk);
277 if (chunk == (unsigned long)-1 && PyErr_Occurred())
278 goto Done;
279 newn = PyNumber_Rshift(n, thirtytwo);
280 if (newn == NULL)
281 goto Done;
282 Py_DECREF(n);
283 n = newn;
284 if (keyused >= keymax) {
285 unsigned long bigger = keymax << 1;
286 if ((bigger >> 1) != keymax) {
287 PyErr_NoMemory();
288 goto Done;
289 }
290 key = (unsigned long *)PyMem_Realloc(key,
291 bigger * sizeof(*key));
292 if (key == NULL)
293 goto Done;
294 keymax = bigger;
295 }
296 assert(keyused < keymax);
297 key[keyused++] = chunk;
298 }
299
300 if (keyused == 0)
301 key[keyused++] = 0UL;
302 result = init_by_array(self, key, keyused);
303 Done:
304 Py_XDECREF(masklower);
305 Py_XDECREF(thirtytwo);
306 Py_XDECREF(n);
307 PyMem_Free(key);
308 return result;
309 }
310
311 static PyObject *
312 random_getstate(RandomObject *self)
313 {
314 PyObject *state;
315 PyObject *element;
316 int i;
317
318 state = PyTuple_New(N+1);
319 if (state == NULL)
320 return NULL;
321 for (i=0; i<N ; i++) {
322 element = PyLong_FromUnsignedLong(self->state[i]);
323 if (element == NULL)
324 goto Fail;
325 PyTuple_SET_ITEM(state, i, element);
326 }
327 element = PyLong_FromLong((long)(self->index));
328 if (element == NULL)
329 goto Fail;
330 PyTuple_SET_ITEM(state, i, element);
331 return state;
332
333 Fail:
334 Py_DECREF(state);
335 return NULL;
336 }
337
338 static PyObject *
339 random_setstate(RandomObject *self, PyObject *state)
340 {
341 int i;
342 unsigned long element;
343 long index;
344
345 if (!PyTuple_Check(state)) {
346 PyErr_SetString(PyExc_TypeError,
347 "state vector must be a tuple");
348 return NULL;
349 }
350 if (PyTuple_Size(state) != N+1) {
351 PyErr_SetString(PyExc_ValueError,
352 "state vector is the wrong size");
353 return NULL;
354 }
355
356 for (i=0; i<N ; i++) {
357 element = PyLong_AsUnsignedLong(PyTuple_GET_ITEM(state, i));
358 if (element == (unsigned long)-1 && PyErr_Occurred())
359 return NULL;
360 self->state[i] = element & 0xffffffffUL; /* Make sure we get sane state */
361 }
362
363 index = PyLong_AsLong(PyTuple_GET_ITEM(state, i));
364 if (index == -1 && PyErr_Occurred())
365 return NULL;
366 self->index = (int)index;
367
368 Py_INCREF(Py_None);
369 return Py_None;
370 }
371
372 /*
373 Jumpahead should be a fast way advance the generator n-steps ahead, but
374 lacking a formula for that, the next best is to use n and the existing
375 state to create a new state far away from the original.
376
377 The generator uses constant spaced additive feedback, so shuffling the
378 state elements ought to produce a state which would not be encountered
379 (in the near term) by calls to random(). Shuffling is normally
380 implemented by swapping the ith element with another element ranging
381 from 0 to i inclusive. That allows the element to have the possibility
382 of not being moved. Since the goal is to produce a new, different
383 state, the swap element is ranged from 0 to i-1 inclusive. This assures
384 that each element gets moved at least once.
385
386 To make sure that consecutive calls to jumpahead(n) produce different
387 states (even in the rare case of involutory shuffles), i+1 is added to
388 each element at position i. Successive calls are then guaranteed to
389 have changing (growing) values as well as shuffled positions.
390
391 Finally, the self->index value is set to N so that the generator itself
392 kicks in on the next call to random(). This assures that all results
393 have been through the generator and do not just reflect alterations to
394 the underlying state.
395 */
396
397 static PyObject *
398 random_jumpahead(RandomObject *self, PyObject *n)
399 {
400 long i, j;
401 PyObject *iobj;
402 PyObject *remobj;
403 unsigned long *mt, tmp;
404
405 if (!PyInt_Check(n) && !PyLong_Check(n)) {
406 PyErr_Format(PyExc_TypeError, "jumpahead requires an "
407 "integer, not '%s'",
408 Py_TYPE(n)->tp_name);
409 return NULL;
410 }
411
412 mt = self->state;
413 for (i = N-1; i > 1; i--) {
414 iobj = PyInt_FromLong(i);
415 if (iobj == NULL)
416 return NULL;
417 remobj = PyNumber_Remainder(n, iobj);
418 Py_DECREF(iobj);
419 if (remobj == NULL)
420 return NULL;
421 j = PyInt_AsLong(remobj);
422 Py_DECREF(remobj);
423 if (j == -1L && PyErr_Occurred())
424 return NULL;
425 tmp = mt[i];
426 mt[i] = mt[j];
427 mt[j] = tmp;
428 }
429
430 for (i = 0; i < N; i++)
431 mt[i] += i+1;
432
433 self->index = N;
434 Py_INCREF(Py_None);
435 return Py_None;
436 }
437
438 static PyObject *
439 random_getrandbits(RandomObject *self, PyObject *args)
440 {
441 int k, i, bytes;
442 unsigned long r;
443 unsigned char *bytearray;
444 PyObject *result;
445
446 if (!PyArg_ParseTuple(args, "i:getrandbits", &k))
447 return NULL;
448
449 if (k <= 0) {
450 PyErr_SetString(PyExc_ValueError,
451 "number of bits must be greater than zero");
452 return NULL;
453 }
454
455 bytes = ((k - 1) / 32 + 1) * 4;
456 bytearray = (unsigned char *)PyMem_Malloc(bytes);
457 if (bytearray == NULL) {
458 PyErr_NoMemory();
459 return NULL;
460 }
461
462 /* Fill-out whole words, byte-by-byte to avoid endianness issues */
463 for (i=0 ; i<bytes ; i+=4, k-=32) {
464 r = genrand_int32(self);
465 if (k < 32)
466 r >>= (32 - k);
467 bytearray[i+0] = (unsigned char)r;
468 bytearray[i+1] = (unsigned char)(r >> 8);
469 bytearray[i+2] = (unsigned char)(r >> 16);
470 bytearray[i+3] = (unsigned char)(r >> 24);
471 }
472
473 /* little endian order to match bytearray assignment order */
474 result = _PyLong_FromByteArray(bytearray, bytes, 1, 0);
475 PyMem_Free(bytearray);
476 return result;
477 }
478
479 static PyObject *
480 random_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
481 {
482 RandomObject *self;
483 PyObject *tmp;
484
485 if (type == &Random_Type && !_PyArg_NoKeywords("Random()", kwds))
486 return NULL;
487
488 self = (RandomObject *)type->tp_alloc(type, 0);
489 if (self == NULL)
490 return NULL;
491 tmp = random_seed(self, args);
492 if (tmp == NULL) {
493 Py_DECREF(self);
494 return NULL;
495 }
496 Py_DECREF(tmp);
497 return (PyObject *)self;
498 }
499
500 static PyMethodDef random_methods[] = {
501 {"random", (PyCFunction)random_random, METH_NOARGS,
502 PyDoc_STR("random() -> x in the interval [0, 1).")},
503 {"seed", (PyCFunction)random_seed, METH_VARARGS,
504 PyDoc_STR("seed([n]) -> None. Defaults to current time.")},
505 {"getstate", (PyCFunction)random_getstate, METH_NOARGS,
506 PyDoc_STR("getstate() -> tuple containing the current state.")},
507 {"setstate", (PyCFunction)random_setstate, METH_O,
508 PyDoc_STR("setstate(state) -> None. Restores generator state.")},
509 {"jumpahead", (PyCFunction)random_jumpahead, METH_O,
510 PyDoc_STR("jumpahead(int) -> None. Create new state from "
511 "existing state and integer.")},
512 {"getrandbits", (PyCFunction)random_getrandbits, METH_VARARGS,
513 PyDoc_STR("getrandbits(k) -> x. Generates a long int with "
514 "k random bits.")},
515 {NULL, NULL} /* sentinel */
516 };
517
518 PyDoc_STRVAR(random_doc,
519 "Random() -> create a random number generator with its own internal state.");
520
521 static PyTypeObject Random_Type = {
522 PyVarObject_HEAD_INIT(NULL, 0)
523 "_random.Random", /*tp_name*/
524 sizeof(RandomObject), /*tp_basicsize*/
525 0, /*tp_itemsize*/
526 /* methods */
527 0, /*tp_dealloc*/
528 0, /*tp_print*/
529 0, /*tp_getattr*/
530 0, /*tp_setattr*/
531 0, /*tp_compare*/
532 0, /*tp_repr*/
533 0, /*tp_as_number*/
534 0, /*tp_as_sequence*/
535 0, /*tp_as_mapping*/
536 0, /*tp_hash*/
537 0, /*tp_call*/
538 0, /*tp_str*/
539 PyObject_GenericGetAttr, /*tp_getattro*/
540 0, /*tp_setattro*/
541 0, /*tp_as_buffer*/
542 Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /*tp_flags*/
543 random_doc, /*tp_doc*/
544 0, /*tp_traverse*/
545 0, /*tp_clear*/
546 0, /*tp_richcompare*/
547 0, /*tp_weaklistoffset*/
548 0, /*tp_iter*/
549 0, /*tp_iternext*/
550 random_methods, /*tp_methods*/
551 0, /*tp_members*/
552 0, /*tp_getset*/
553 0, /*tp_base*/
554 0, /*tp_dict*/
555 0, /*tp_descr_get*/
556 0, /*tp_descr_set*/
557 0, /*tp_dictoffset*/
558 0, /*tp_init*/
559 0, /*tp_alloc*/
560 random_new, /*tp_new*/
561 _PyObject_Del, /*tp_free*/
562 0, /*tp_is_gc*/
563 };
564
565 PyDoc_STRVAR(module_doc,
566 "Module implements the Mersenne Twister random number generator.");
567
568 PyMODINIT_FUNC
569 init_random(void)
570 {
571 PyObject *m;
572
573 if (PyType_Ready(&Random_Type) < 0)
574 return;
575 m = Py_InitModule3("_random", NULL, module_doc);
576 if (m == NULL)
577 return;
578 Py_INCREF(&Random_Type);
579 PyModule_AddObject(m, "Random", (PyObject *)&Random_Type);
580 }