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