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_uint32() 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 "pycore_byteswap.h" // _Py_bswap32()
71 #ifdef HAVE_PROCESS_H
72 # include <process.h> // 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 *Random_Type;
84 PyObject *Long___abs__;
85 } _randomstate;
86
87 static inline _randomstate*
get_random_state(PyObject * module)88 get_random_state(PyObject *module)
89 {
90 void *state = PyModule_GetState(module);
91 assert(state != NULL);
92 return (_randomstate *)state;
93 }
94
95 static struct PyModuleDef _randommodule;
96
97 #define _randomstate_global get_random_state(PyState_FindModule(&_randommodule))
98
99 typedef struct {
100 PyObject_HEAD
101 int index;
102 uint32_t state[N];
103 } RandomObject;
104
105
106 #include "clinic/_randommodule.c.h"
107
108 /*[clinic input]
109 module _random
110 class _random.Random "RandomObject *" "&Random_Type"
111 [clinic start generated code]*/
112 /*[clinic end generated code: output=da39a3ee5e6b4b0d input=f79898ae7847c321]*/
113
114 /* Random methods */
115
116
117 /* generates a random number on [0,0xffffffff]-interval */
118 static uint32_t
genrand_uint32(RandomObject * self)119 genrand_uint32(RandomObject *self)
120 {
121 uint32_t y;
122 static const uint32_t mag01[2] = {0x0U, MATRIX_A};
123 /* mag01[x] = x * MATRIX_A for x=0,1 */
124 uint32_t *mt;
125
126 mt = self->state;
127 if (self->index >= N) { /* generate N words at one time */
128 int kk;
129
130 for (kk=0;kk<N-M;kk++) {
131 y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
132 mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1U];
133 }
134 for (;kk<N-1;kk++) {
135 y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
136 mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1U];
137 }
138 y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
139 mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1U];
140
141 self->index = 0;
142 }
143
144 y = mt[self->index++];
145 y ^= (y >> 11);
146 y ^= (y << 7) & 0x9d2c5680U;
147 y ^= (y << 15) & 0xefc60000U;
148 y ^= (y >> 18);
149 return y;
150 }
151
152 /* random_random is the function named genrand_res53 in the original code;
153 * generates a random number on [0,1) with 53-bit resolution; note that
154 * 9007199254740992 == 2**53; I assume they're spelling "/2**53" as
155 * multiply-by-reciprocal in the (likely vain) hope that the compiler will
156 * optimize the division away at compile-time. 67108864 is 2**26. In
157 * effect, a contains 27 random bits shifted left 26, and b fills in the
158 * lower 26 bits of the 53-bit numerator.
159 * The original code credited Isaku Wada for this algorithm, 2002/01/09.
160 */
161
162 /*[clinic input]
163 _random.Random.random
164
165 self: self(type="RandomObject *")
166
167 random() -> x in the interval [0, 1).
168 [clinic start generated code]*/
169
170 static PyObject *
_random_Random_random_impl(RandomObject * self)171 _random_Random_random_impl(RandomObject *self)
172 /*[clinic end generated code: output=117ff99ee53d755c input=afb2a59cbbb00349]*/
173 {
174 uint32_t a=genrand_uint32(self)>>5, b=genrand_uint32(self)>>6;
175 return PyFloat_FromDouble((a*67108864.0+b)*(1.0/9007199254740992.0));
176 }
177
178 /* initializes mt[N] with a seed */
179 static void
init_genrand(RandomObject * self,uint32_t s)180 init_genrand(RandomObject *self, uint32_t s)
181 {
182 int mti;
183 uint32_t *mt;
184
185 mt = self->state;
186 mt[0]= s;
187 for (mti=1; mti<N; mti++) {
188 mt[mti] =
189 (1812433253U * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
190 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
191 /* In the previous versions, MSBs of the seed affect */
192 /* only MSBs of the array mt[]. */
193 /* 2002/01/09 modified by Makoto Matsumoto */
194 }
195 self->index = mti;
196 return;
197 }
198
199 /* initialize by an array with array-length */
200 /* init_key is the array for initializing keys */
201 /* key_length is its length */
202 static void
init_by_array(RandomObject * self,uint32_t init_key[],size_t key_length)203 init_by_array(RandomObject *self, uint32_t init_key[], size_t key_length)
204 {
205 size_t i, j, k; /* was signed in the original code. RDH 12/16/2002 */
206 uint32_t *mt;
207
208 mt = self->state;
209 init_genrand(self, 19650218U);
210 i=1; j=0;
211 k = (N>key_length ? N : key_length);
212 for (; k; k--) {
213 mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525U))
214 + init_key[j] + (uint32_t)j; /* non linear */
215 i++; j++;
216 if (i>=N) { mt[0] = mt[N-1]; i=1; }
217 if (j>=key_length) j=0;
218 }
219 for (k=N-1; k; k--) {
220 mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941U))
221 - (uint32_t)i; /* non linear */
222 i++;
223 if (i>=N) { mt[0] = mt[N-1]; i=1; }
224 }
225
226 mt[0] = 0x80000000U; /* MSB is 1; assuring non-zero initial array */
227 }
228
229 /*
230 * The rest is Python-specific code, neither part of, nor derived from, the
231 * Twister download.
232 */
233
234 static int
random_seed_urandom(RandomObject * self)235 random_seed_urandom(RandomObject *self)
236 {
237 uint32_t key[N];
238
239 if (_PyOS_URandomNonblock(key, sizeof(key)) < 0) {
240 return -1;
241 }
242 init_by_array(self, key, Py_ARRAY_LENGTH(key));
243 return 0;
244 }
245
246 static void
random_seed_time_pid(RandomObject * self)247 random_seed_time_pid(RandomObject *self)
248 {
249 _PyTime_t now;
250 uint32_t key[5];
251
252 now = _PyTime_GetSystemClock();
253 key[0] = (uint32_t)(now & 0xffffffffU);
254 key[1] = (uint32_t)(now >> 32);
255
256 key[2] = (uint32_t)getpid();
257
258 now = _PyTime_GetMonotonicClock();
259 key[3] = (uint32_t)(now & 0xffffffffU);
260 key[4] = (uint32_t)(now >> 32);
261
262 init_by_array(self, key, Py_ARRAY_LENGTH(key));
263 }
264
265 static PyObject *
random_seed(RandomObject * self,PyObject * arg)266 random_seed(RandomObject *self, PyObject *arg)
267 {
268 PyObject *result = NULL; /* guilty until proved innocent */
269 PyObject *n = NULL;
270 uint32_t *key = NULL;
271 size_t bits, keyused;
272 int res;
273
274 if (arg == NULL || arg == Py_None) {
275 if (random_seed_urandom(self) < 0) {
276 PyErr_Clear();
277
278 /* Reading system entropy failed, fall back on the worst entropy:
279 use the current time and process identifier. */
280 random_seed_time_pid(self);
281 }
282 Py_RETURN_NONE;
283 }
284
285 /* This algorithm relies on the number being unsigned.
286 * So: if the arg is a PyLong, use its absolute value.
287 * Otherwise use its hash value, cast to unsigned.
288 */
289 if (PyLong_CheckExact(arg)) {
290 n = PyNumber_Absolute(arg);
291 } else if (PyLong_Check(arg)) {
292 /* Calling int.__abs__() prevents calling arg.__abs__(), which might
293 return an invalid value. See issue #31478. */
294 n = PyObject_CallOneArg(_randomstate_global->Long___abs__, arg);
295 }
296 else {
297 Py_hash_t hash = PyObject_Hash(arg);
298 if (hash == -1)
299 goto Done;
300 n = PyLong_FromSize_t((size_t)hash);
301 }
302 if (n == NULL)
303 goto Done;
304
305 /* Now split n into 32-bit chunks, from the right. */
306 bits = _PyLong_NumBits(n);
307 if (bits == (size_t)-1 && PyErr_Occurred())
308 goto Done;
309
310 /* Figure out how many 32-bit chunks this gives us. */
311 keyused = bits == 0 ? 1 : (bits - 1) / 32 + 1;
312
313 /* Convert seed to byte sequence. */
314 key = (uint32_t *)PyMem_Malloc((size_t)4 * keyused);
315 if (key == NULL) {
316 PyErr_NoMemory();
317 goto Done;
318 }
319 res = _PyLong_AsByteArray((PyLongObject *)n,
320 (unsigned char *)key, keyused * 4,
321 PY_LITTLE_ENDIAN,
322 0); /* unsigned */
323 if (res == -1) {
324 goto Done;
325 }
326
327 #if PY_BIG_ENDIAN
328 {
329 size_t i, j;
330 /* Reverse an array. */
331 for (i = 0, j = keyused - 1; i < j; i++, j--) {
332 uint32_t tmp = key[i];
333 key[i] = key[j];
334 key[j] = tmp;
335 }
336 }
337 #endif
338 init_by_array(self, key, keyused);
339
340 Py_INCREF(Py_None);
341 result = Py_None;
342
343 Done:
344 Py_XDECREF(n);
345 PyMem_Free(key);
346 return result;
347 }
348
349 /*[clinic input]
350 _random.Random.seed
351
352 self: self(type="RandomObject *")
353 n: object = None
354 /
355
356 seed([n]) -> None.
357
358 Defaults to use urandom and falls back to a combination
359 of the current time and the process identifier.
360 [clinic start generated code]*/
361
362 static PyObject *
_random_Random_seed_impl(RandomObject * self,PyObject * n)363 _random_Random_seed_impl(RandomObject *self, PyObject *n)
364 /*[clinic end generated code: output=0fad1e16ba883681 input=78d6ef0d52532a54]*/
365 {
366 return random_seed(self, n);
367 }
368
369 /*[clinic input]
370 _random.Random.getstate
371
372 self: self(type="RandomObject *")
373
374 getstate() -> tuple containing the current state.
375 [clinic start generated code]*/
376
377 static PyObject *
_random_Random_getstate_impl(RandomObject * self)378 _random_Random_getstate_impl(RandomObject *self)
379 /*[clinic end generated code: output=bf6cef0c092c7180 input=b937a487928c0e89]*/
380 {
381 PyObject *state;
382 PyObject *element;
383 int i;
384
385 state = PyTuple_New(N+1);
386 if (state == NULL)
387 return NULL;
388 for (i=0; i<N ; i++) {
389 element = PyLong_FromUnsignedLong(self->state[i]);
390 if (element == NULL)
391 goto Fail;
392 PyTuple_SET_ITEM(state, i, element);
393 }
394 element = PyLong_FromLong((long)(self->index));
395 if (element == NULL)
396 goto Fail;
397 PyTuple_SET_ITEM(state, i, element);
398 return state;
399
400 Fail:
401 Py_DECREF(state);
402 return NULL;
403 }
404
405
406 /*[clinic input]
407 _random.Random.setstate
408
409 self: self(type="RandomObject *")
410 state: object
411 /
412
413 setstate(state) -> None. Restores generator state.
414 [clinic start generated code]*/
415
416 static PyObject *
_random_Random_setstate(RandomObject * self,PyObject * state)417 _random_Random_setstate(RandomObject *self, PyObject *state)
418 /*[clinic end generated code: output=fd1c3cd0037b6681 input=b3b4efbb1bc66af8]*/
419 {
420 int i;
421 unsigned long element;
422 long index;
423 uint32_t new_state[N];
424
425 if (!PyTuple_Check(state)) {
426 PyErr_SetString(PyExc_TypeError,
427 "state vector must be a tuple");
428 return NULL;
429 }
430 if (PyTuple_Size(state) != N+1) {
431 PyErr_SetString(PyExc_ValueError,
432 "state vector is the wrong size");
433 return NULL;
434 }
435
436 for (i=0; i<N ; i++) {
437 element = PyLong_AsUnsignedLong(PyTuple_GET_ITEM(state, i));
438 if (element == (unsigned long)-1 && PyErr_Occurred())
439 return NULL;
440 new_state[i] = (uint32_t)element;
441 }
442
443 index = PyLong_AsLong(PyTuple_GET_ITEM(state, i));
444 if (index == -1 && PyErr_Occurred())
445 return NULL;
446 if (index < 0 || index > N) {
447 PyErr_SetString(PyExc_ValueError, "invalid state");
448 return NULL;
449 }
450 self->index = (int)index;
451 for (i = 0; i < N; i++)
452 self->state[i] = new_state[i];
453
454 Py_RETURN_NONE;
455 }
456
457 /*[clinic input]
458
459 _random.Random.getrandbits
460
461 self: self(type="RandomObject *")
462 k: int
463 /
464
465 getrandbits(k) -> x. Generates an int with k random bits.
466 [clinic start generated code]*/
467
468 static PyObject *
_random_Random_getrandbits_impl(RandomObject * self,int k)469 _random_Random_getrandbits_impl(RandomObject *self, int k)
470 /*[clinic end generated code: output=b402f82a2158887f input=8c0e6396dd176fc0]*/
471 {
472 int i, words;
473 uint32_t r;
474 uint32_t *wordarray;
475 PyObject *result;
476
477 if (k < 0) {
478 PyErr_SetString(PyExc_ValueError,
479 "number of bits must be non-negative");
480 return NULL;
481 }
482
483 if (k == 0)
484 return PyLong_FromLong(0);
485
486 if (k <= 32) /* Fast path */
487 return PyLong_FromUnsignedLong(genrand_uint32(self) >> (32 - k));
488
489 words = (k - 1) / 32 + 1;
490 wordarray = (uint32_t *)PyMem_Malloc(words * 4);
491 if (wordarray == NULL) {
492 PyErr_NoMemory();
493 return NULL;
494 }
495
496 /* Fill-out bits of long integer, by 32-bit words, from least significant
497 to most significant. */
498 #if PY_LITTLE_ENDIAN
499 for (i = 0; i < words; i++, k -= 32)
500 #else
501 for (i = words - 1; i >= 0; i--, k -= 32)
502 #endif
503 {
504 r = genrand_uint32(self);
505 if (k < 32)
506 r >>= (32 - k); /* Drop least significant bits */
507 wordarray[i] = r;
508 }
509
510 result = _PyLong_FromByteArray((unsigned char *)wordarray, words * 4,
511 PY_LITTLE_ENDIAN, 0 /* unsigned */);
512 PyMem_Free(wordarray);
513 return result;
514 }
515
516 static PyObject *
random_new(PyTypeObject * type,PyObject * args,PyObject * kwds)517 random_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
518 {
519 RandomObject *self;
520 PyObject *tmp;
521
522 if (type == (PyTypeObject*)_randomstate_global->Random_Type &&
523 !_PyArg_NoKeywords("Random()", kwds)) {
524 return NULL;
525 }
526
527 self = (RandomObject *)PyType_GenericAlloc(type, 0);
528 if (self == NULL)
529 return NULL;
530 tmp = random_seed(self, args);
531 if (tmp == NULL) {
532 Py_DECREF(self);
533 return NULL;
534 }
535 Py_DECREF(tmp);
536 return (PyObject *)self;
537 }
538
539
540 static PyMethodDef random_methods[] = {
541 _RANDOM_RANDOM_RANDOM_METHODDEF
542 _RANDOM_RANDOM_SEED_METHODDEF
543 _RANDOM_RANDOM_GETSTATE_METHODDEF
544 _RANDOM_RANDOM_SETSTATE_METHODDEF
545 _RANDOM_RANDOM_GETRANDBITS_METHODDEF
546 {NULL, NULL} /* sentinel */
547 };
548
549 PyDoc_STRVAR(random_doc,
550 "Random() -> create a random number generator with its own internal state.");
551
552 static PyType_Slot Random_Type_slots[] = {
553 {Py_tp_doc, (void *)random_doc},
554 {Py_tp_methods, random_methods},
555 {Py_tp_new, random_new},
556 {Py_tp_free, PyObject_Free},
557 {0, 0},
558 };
559
560 static PyType_Spec Random_Type_spec = {
561 "_random.Random",
562 sizeof(RandomObject),
563 0,
564 Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE,
565 Random_Type_slots
566 };
567
568 PyDoc_STRVAR(module_doc,
569 "Module implements the Mersenne Twister random number generator.");
570
571 static int
_random_traverse(PyObject * module,visitproc visit,void * arg)572 _random_traverse(PyObject *module, visitproc visit, void *arg)
573 {
574 Py_VISIT(get_random_state(module)->Random_Type);
575 return 0;
576 }
577
578 static int
_random_clear(PyObject * module)579 _random_clear(PyObject *module)
580 {
581 Py_CLEAR(get_random_state(module)->Random_Type);
582 Py_CLEAR(get_random_state(module)->Long___abs__);
583 return 0;
584 }
585
586 static void
_random_free(void * module)587 _random_free(void *module)
588 {
589 _random_clear((PyObject *)module);
590 }
591
592 static struct PyModuleDef _randommodule = {
593 PyModuleDef_HEAD_INIT,
594 "_random",
595 module_doc,
596 sizeof(_randomstate),
597 NULL,
598 NULL,
599 _random_traverse,
600 _random_clear,
601 _random_free,
602 };
603
604 PyMODINIT_FUNC
PyInit__random(void)605 PyInit__random(void)
606 {
607 PyObject *m;
608
609 PyObject *Random_Type = PyType_FromSpec(&Random_Type_spec);
610 if (Random_Type == NULL) {
611 return NULL;
612 }
613
614 m = PyModule_Create(&_randommodule);
615 if (m == NULL) {
616 Py_DECREF(Random_Type);
617 return NULL;
618 }
619 get_random_state(m)->Random_Type = Random_Type;
620
621 Py_INCREF(Random_Type);
622 PyModule_AddObject(m, "Random", Random_Type);
623
624 /* Look up and save int.__abs__, which is needed in random_seed(). */
625 PyObject *longval = NULL, *longtype = NULL;
626 longval = PyLong_FromLong(0);
627 if (longval == NULL) goto fail;
628
629 longtype = PyObject_Type(longval);
630 if (longtype == NULL) goto fail;
631
632 PyObject *abs = PyObject_GetAttrString(longtype, "__abs__");
633 if (abs == NULL) goto fail;
634
635 Py_DECREF(longtype);
636 Py_DECREF(longval);
637 get_random_state(m)->Long___abs__ = abs;
638
639 return m;
640
641 fail:
642 Py_XDECREF(longtype);
643 Py_XDECREF(longval);
644 Py_DECREF(m);
645 return NULL;
646 }
647