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 #ifdef HAVE_PROCESS_H
72 # include <process.h> /* needed for getpid() */
73 #endif
74
75 /* Period parameters -- These are all magic. Don't change. */
76 #define N 624
77 #define M 397
78 #define MATRIX_A 0x9908b0dfU /* constant vector a */
79 #define UPPER_MASK 0x80000000U /* most significant w-r bits */
80 #define LOWER_MASK 0x7fffffffU /* least significant r bits */
81
82 typedef struct {
83 PyObject_HEAD
84 int index;
85 uint32_t state[N];
86 } RandomObject;
87
88 static PyTypeObject Random_Type;
89
90 #define RandomObject_Check(v) (Py_TYPE(v) == &Random_Type)
91
92
93 /* Random methods */
94
95
96 /* generates a random number on [0,0xffffffff]-interval */
97 static uint32_t
genrand_int32(RandomObject * self)98 genrand_int32(RandomObject *self)
99 {
100 uint32_t y;
101 static const uint32_t mag01[2] = {0x0U, MATRIX_A};
102 /* mag01[x] = x * MATRIX_A for x=0,1 */
103 uint32_t *mt;
104
105 mt = self->state;
106 if (self->index >= N) { /* generate N words at one time */
107 int kk;
108
109 for (kk=0;kk<N-M;kk++) {
110 y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
111 mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1U];
112 }
113 for (;kk<N-1;kk++) {
114 y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
115 mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1U];
116 }
117 y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
118 mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1U];
119
120 self->index = 0;
121 }
122
123 y = mt[self->index++];
124 y ^= (y >> 11);
125 y ^= (y << 7) & 0x9d2c5680U;
126 y ^= (y << 15) & 0xefc60000U;
127 y ^= (y >> 18);
128 return y;
129 }
130
131 /* random_random is the function named genrand_res53 in the original code;
132 * generates a random number on [0,1) with 53-bit resolution; note that
133 * 9007199254740992 == 2**53; I assume they're spelling "/2**53" as
134 * multiply-by-reciprocal in the (likely vain) hope that the compiler will
135 * optimize the division away at compile-time. 67108864 is 2**26. In
136 * effect, a contains 27 random bits shifted left 26, and b fills in the
137 * lower 26 bits of the 53-bit numerator.
138 * The original code credited Isaku Wada for this algorithm, 2002/01/09.
139 */
140 static PyObject *
random_random(RandomObject * self)141 random_random(RandomObject *self)
142 {
143 uint32_t a=genrand_int32(self)>>5, b=genrand_int32(self)>>6;
144 return PyFloat_FromDouble((a*67108864.0+b)*(1.0/9007199254740992.0));
145 }
146
147 /* initializes mt[N] with a seed */
148 static void
init_genrand(RandomObject * self,uint32_t s)149 init_genrand(RandomObject *self, uint32_t s)
150 {
151 int mti;
152 uint32_t *mt;
153
154 mt = self->state;
155 mt[0]= s;
156 for (mti=1; mti<N; mti++) {
157 mt[mti] =
158 (1812433253U * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
159 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
160 /* In the previous versions, MSBs of the seed affect */
161 /* only MSBs of the array mt[]. */
162 /* 2002/01/09 modified by Makoto Matsumoto */
163 }
164 self->index = mti;
165 return;
166 }
167
168 /* initialize by an array with array-length */
169 /* init_key is the array for initializing keys */
170 /* key_length is its length */
171 static void
init_by_array(RandomObject * self,uint32_t init_key[],size_t key_length)172 init_by_array(RandomObject *self, uint32_t init_key[], size_t key_length)
173 {
174 size_t i, j, k; /* was signed in the original code. RDH 12/16/2002 */
175 uint32_t *mt;
176
177 mt = self->state;
178 init_genrand(self, 19650218U);
179 i=1; j=0;
180 k = (N>key_length ? N : key_length);
181 for (; k; k--) {
182 mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525U))
183 + init_key[j] + (uint32_t)j; /* non linear */
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)) * 1566083941U))
190 - (uint32_t)i; /* non linear */
191 i++;
192 if (i>=N) { mt[0] = mt[N-1]; i=1; }
193 }
194
195 mt[0] = 0x80000000U; /* MSB is 1; assuring non-zero initial array */
196 }
197
198 /*
199 * The rest is Python-specific code, neither part of, nor derived from, the
200 * Twister download.
201 */
202
203 static int
random_seed_urandom(RandomObject * self)204 random_seed_urandom(RandomObject *self)
205 {
206 PY_UINT32_T key[N];
207
208 if (_PyOS_URandomNonblock(key, sizeof(key)) < 0) {
209 return -1;
210 }
211 init_by_array(self, key, Py_ARRAY_LENGTH(key));
212 return 0;
213 }
214
215 static void
random_seed_time_pid(RandomObject * self)216 random_seed_time_pid(RandomObject *self)
217 {
218 _PyTime_t now;
219 uint32_t key[5];
220
221 now = _PyTime_GetSystemClock();
222 key[0] = (PY_UINT32_T)(now & 0xffffffffU);
223 key[1] = (PY_UINT32_T)(now >> 32);
224
225 key[2] = (PY_UINT32_T)getpid();
226
227 now = _PyTime_GetMonotonicClock();
228 key[3] = (PY_UINT32_T)(now & 0xffffffffU);
229 key[4] = (PY_UINT32_T)(now >> 32);
230
231 init_by_array(self, key, Py_ARRAY_LENGTH(key));
232 }
233
234 static PyObject *
random_seed(RandomObject * self,PyObject * args)235 random_seed(RandomObject *self, PyObject *args)
236 {
237 PyObject *result = NULL; /* guilty until proved innocent */
238 PyObject *n = NULL;
239 uint32_t *key = NULL;
240 size_t bits, keyused;
241 int res;
242 PyObject *arg = NULL;
243
244 if (!PyArg_UnpackTuple(args, "seed", 0, 1, &arg))
245 return NULL;
246
247 if (arg == NULL || arg == Py_None) {
248 if (random_seed_urandom(self) < 0) {
249 PyErr_Clear();
250
251 /* Reading system entropy failed, fall back on the worst entropy:
252 use the current time and process identifier. */
253 random_seed_time_pid(self);
254 }
255 Py_RETURN_NONE;
256 }
257
258 /* This algorithm relies on the number being unsigned.
259 * So: if the arg is a PyLong, use its absolute value.
260 * Otherwise use its hash value, cast to unsigned.
261 */
262 if (PyLong_Check(arg)) {
263 /* Calling int.__abs__() prevents calling arg.__abs__(), which might
264 return an invalid value. See issue #31478. */
265 n = PyLong_Type.tp_as_number->nb_absolute(arg);
266 }
267 else {
268 Py_hash_t hash = PyObject_Hash(arg);
269 if (hash == -1)
270 goto Done;
271 n = PyLong_FromSize_t((size_t)hash);
272 }
273 if (n == NULL)
274 goto Done;
275
276 /* Now split n into 32-bit chunks, from the right. */
277 bits = _PyLong_NumBits(n);
278 if (bits == (size_t)-1 && PyErr_Occurred())
279 goto Done;
280
281 /* Figure out how many 32-bit chunks this gives us. */
282 keyused = bits == 0 ? 1 : (bits - 1) / 32 + 1;
283
284 /* Convert seed to byte sequence. */
285 key = (uint32_t *)PyMem_Malloc((size_t)4 * keyused);
286 if (key == NULL) {
287 PyErr_NoMemory();
288 goto Done;
289 }
290 res = _PyLong_AsByteArray((PyLongObject *)n,
291 (unsigned char *)key, keyused * 4,
292 PY_LITTLE_ENDIAN,
293 0); /* unsigned */
294 if (res == -1) {
295 goto Done;
296 }
297
298 #if PY_BIG_ENDIAN
299 {
300 size_t i, j;
301 /* Reverse an array. */
302 for (i = 0, j = keyused - 1; i < j; i++, j--) {
303 uint32_t tmp = key[i];
304 key[i] = key[j];
305 key[j] = tmp;
306 }
307 }
308 #endif
309 init_by_array(self, key, keyused);
310
311 Py_INCREF(Py_None);
312 result = Py_None;
313
314 Done:
315 Py_XDECREF(n);
316 PyMem_Free(key);
317 return result;
318 }
319
320 static PyObject *
random_getstate(RandomObject * self)321 random_getstate(RandomObject *self)
322 {
323 PyObject *state;
324 PyObject *element;
325 int i;
326
327 state = PyTuple_New(N+1);
328 if (state == NULL)
329 return NULL;
330 for (i=0; i<N ; i++) {
331 element = PyLong_FromUnsignedLong(self->state[i]);
332 if (element == NULL)
333 goto Fail;
334 PyTuple_SET_ITEM(state, i, element);
335 }
336 element = PyLong_FromLong((long)(self->index));
337 if (element == NULL)
338 goto Fail;
339 PyTuple_SET_ITEM(state, i, element);
340 return state;
341
342 Fail:
343 Py_DECREF(state);
344 return NULL;
345 }
346
347 static PyObject *
random_setstate(RandomObject * self,PyObject * state)348 random_setstate(RandomObject *self, PyObject *state)
349 {
350 int i;
351 unsigned long element;
352 long index;
353 uint32_t new_state[N];
354
355 if (!PyTuple_Check(state)) {
356 PyErr_SetString(PyExc_TypeError,
357 "state vector must be a tuple");
358 return NULL;
359 }
360 if (PyTuple_Size(state) != N+1) {
361 PyErr_SetString(PyExc_ValueError,
362 "state vector is the wrong size");
363 return NULL;
364 }
365
366 for (i=0; i<N ; i++) {
367 element = PyLong_AsUnsignedLong(PyTuple_GET_ITEM(state, i));
368 if (element == (unsigned long)-1 && PyErr_Occurred())
369 return NULL;
370 new_state[i] = (uint32_t)element;
371 }
372
373 index = PyLong_AsLong(PyTuple_GET_ITEM(state, i));
374 if (index == -1 && PyErr_Occurred())
375 return NULL;
376 if (index < 0 || index > N) {
377 PyErr_SetString(PyExc_ValueError, "invalid state");
378 return NULL;
379 }
380 self->index = (int)index;
381 for (i = 0; i < N; i++)
382 self->state[i] = new_state[i];
383
384 Py_RETURN_NONE;
385 }
386
387 static PyObject *
random_getrandbits(RandomObject * self,PyObject * args)388 random_getrandbits(RandomObject *self, PyObject *args)
389 {
390 int k, i, words;
391 uint32_t r;
392 uint32_t *wordarray;
393 PyObject *result;
394
395 if (!PyArg_ParseTuple(args, "i:getrandbits", &k))
396 return NULL;
397
398 if (k <= 0) {
399 PyErr_SetString(PyExc_ValueError,
400 "number of bits must be greater than zero");
401 return NULL;
402 }
403
404 if (k <= 32) /* Fast path */
405 return PyLong_FromUnsignedLong(genrand_int32(self) >> (32 - k));
406
407 words = (k - 1) / 32 + 1;
408 wordarray = (uint32_t *)PyMem_Malloc(words * 4);
409 if (wordarray == NULL) {
410 PyErr_NoMemory();
411 return NULL;
412 }
413
414 /* Fill-out bits of long integer, by 32-bit words, from least significant
415 to most significant. */
416 #if PY_LITTLE_ENDIAN
417 for (i = 0; i < words; i++, k -= 32)
418 #else
419 for (i = words - 1; i >= 0; i--, k -= 32)
420 #endif
421 {
422 r = genrand_int32(self);
423 if (k < 32)
424 r >>= (32 - k); /* Drop least significant bits */
425 wordarray[i] = r;
426 }
427
428 result = _PyLong_FromByteArray((unsigned char *)wordarray, words * 4,
429 PY_LITTLE_ENDIAN, 0 /* unsigned */);
430 PyMem_Free(wordarray);
431 return result;
432 }
433
434 static PyObject *
random_new(PyTypeObject * type,PyObject * args,PyObject * kwds)435 random_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
436 {
437 RandomObject *self;
438 PyObject *tmp;
439
440 if (type == &Random_Type && !_PyArg_NoKeywords("Random", kwds))
441 return NULL;
442
443 self = (RandomObject *)type->tp_alloc(type, 0);
444 if (self == NULL)
445 return NULL;
446 tmp = random_seed(self, args);
447 if (tmp == NULL) {
448 Py_DECREF(self);
449 return NULL;
450 }
451 Py_DECREF(tmp);
452 return (PyObject *)self;
453 }
454
455 static PyMethodDef random_methods[] = {
456 {"random", (PyCFunction)random_random, METH_NOARGS,
457 PyDoc_STR("random() -> x in the interval [0, 1).")},
458 {"seed", (PyCFunction)random_seed, METH_VARARGS,
459 PyDoc_STR("seed([n]) -> None. Defaults to current time.")},
460 {"getstate", (PyCFunction)random_getstate, METH_NOARGS,
461 PyDoc_STR("getstate() -> tuple containing the current state.")},
462 {"setstate", (PyCFunction)random_setstate, METH_O,
463 PyDoc_STR("setstate(state) -> None. Restores generator state.")},
464 {"getrandbits", (PyCFunction)random_getrandbits, METH_VARARGS,
465 PyDoc_STR("getrandbits(k) -> x. Generates an int with "
466 "k random bits.")},
467 {NULL, NULL} /* sentinel */
468 };
469
470 PyDoc_STRVAR(random_doc,
471 "Random() -> create a random number generator with its own internal state.");
472
473 static PyTypeObject Random_Type = {
474 PyVarObject_HEAD_INIT(NULL, 0)
475 "_random.Random", /*tp_name*/
476 sizeof(RandomObject), /*tp_basicsize*/
477 0, /*tp_itemsize*/
478 /* methods */
479 0, /*tp_dealloc*/
480 0, /*tp_print*/
481 0, /*tp_getattr*/
482 0, /*tp_setattr*/
483 0, /*tp_reserved*/
484 0, /*tp_repr*/
485 0, /*tp_as_number*/
486 0, /*tp_as_sequence*/
487 0, /*tp_as_mapping*/
488 0, /*tp_hash*/
489 0, /*tp_call*/
490 0, /*tp_str*/
491 PyObject_GenericGetAttr, /*tp_getattro*/
492 0, /*tp_setattro*/
493 0, /*tp_as_buffer*/
494 Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /*tp_flags*/
495 random_doc, /*tp_doc*/
496 0, /*tp_traverse*/
497 0, /*tp_clear*/
498 0, /*tp_richcompare*/
499 0, /*tp_weaklistoffset*/
500 0, /*tp_iter*/
501 0, /*tp_iternext*/
502 random_methods, /*tp_methods*/
503 0, /*tp_members*/
504 0, /*tp_getset*/
505 0, /*tp_base*/
506 0, /*tp_dict*/
507 0, /*tp_descr_get*/
508 0, /*tp_descr_set*/
509 0, /*tp_dictoffset*/
510 0, /*tp_init*/
511 0, /*tp_alloc*/
512 random_new, /*tp_new*/
513 PyObject_Free, /*tp_free*/
514 0, /*tp_is_gc*/
515 };
516
517 PyDoc_STRVAR(module_doc,
518 "Module implements the Mersenne Twister random number generator.");
519
520
521 static struct PyModuleDef _randommodule = {
522 PyModuleDef_HEAD_INIT,
523 "_random",
524 module_doc,
525 -1,
526 NULL,
527 NULL,
528 NULL,
529 NULL,
530 NULL
531 };
532
533 PyMODINIT_FUNC
PyInit__random(void)534 PyInit__random(void)
535 {
536 PyObject *m;
537
538 if (PyType_Ready(&Random_Type) < 0)
539 return NULL;
540 m = PyModule_Create(&_randommodule);
541 if (m == NULL)
542 return NULL;
543 Py_INCREF(&Random_Type);
544 PyModule_AddObject(m, "Random", (PyObject *)&Random_Type);
545 return m;
546 }
547