ยปCore Development>Code coverage>Modules/_randommodule.c

Python code coverage for Modules/_randommodule.c

#countcontent
1n/a/* Random objects */
2n/a
3n/a/* ------------------------------------------------------------------
4n/a The code in this module was based on a download from:
5n/a http://www.math.keio.ac.jp/~matumoto/MT2002/emt19937ar.html
6n/a
7n/a It was modified in 2002 by Raymond Hettinger as follows:
8n/a
9n/a * the principal computational lines untouched.
10n/a
11n/a * renamed genrand_res53() to random_random() and wrapped
12n/a in python calling/return code.
13n/a
14n/a * genrand_int32() and the helper functions, init_genrand()
15n/a and init_by_array(), were declared static, wrapped in
16n/a Python calling/return code. also, their global data
17n/a references were replaced with structure references.
18n/a
19n/a * unused functions from the original were deleted.
20n/a new, original C python code was added to implement the
21n/a Random() interface.
22n/a
23n/a The following are the verbatim comments from the original code:
24n/a
25n/a A C-program for MT19937, with initialization improved 2002/1/26.
26n/a Coded by Takuji Nishimura and Makoto Matsumoto.
27n/a
28n/a Before using, initialize the state by using init_genrand(seed)
29n/a or init_by_array(init_key, key_length).
30n/a
31n/a Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
32n/a All rights reserved.
33n/a
34n/a Redistribution and use in source and binary forms, with or without
35n/a modification, are permitted provided that the following conditions
36n/a are met:
37n/a
38n/a 1. Redistributions of source code must retain the above copyright
39n/a notice, this list of conditions and the following disclaimer.
40n/a
41n/a 2. Redistributions in binary form must reproduce the above copyright
42n/a notice, this list of conditions and the following disclaimer in the
43n/a documentation and/or other materials provided with the distribution.
44n/a
45n/a 3. The names of its contributors may not be used to endorse or promote
46n/a products derived from this software without specific prior written
47n/a permission.
48n/a
49n/a THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
50n/a "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
51n/a LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
52n/a A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
53n/a CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
54n/a EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
55n/a PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
56n/a PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
57n/a LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
58n/a NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
59n/a SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
60n/a
61n/a
62n/a Any feedback is very welcome.
63n/a http://www.math.keio.ac.jp/matumoto/emt.html
64n/a email: matumoto@math.keio.ac.jp
65n/a*/
66n/a
67n/a/* ---------------------------------------------------------------*/
68n/a
69n/a#include "Python.h"
70n/a#include <time.h> /* for seeding to current time */
71n/a#ifdef HAVE_PROCESS_H
72n/a# include <process.h> /* needed for getpid() */
73n/a#endif
74n/a
75n/a/* Period parameters -- These are all magic. Don't change. */
76n/a#define N 624
77n/a#define M 397
78n/a#define MATRIX_A 0x9908b0dfU /* constant vector a */
79n/a#define UPPER_MASK 0x80000000U /* most significant w-r bits */
80n/a#define LOWER_MASK 0x7fffffffU /* least significant r bits */
81n/a
82n/atypedef struct {
83n/a PyObject_HEAD
84n/a int index;
85n/a uint32_t state[N];
86n/a} RandomObject;
87n/a
88n/astatic PyTypeObject Random_Type;
89n/a
90n/a#define RandomObject_Check(v) (Py_TYPE(v) == &Random_Type)
91n/a
92n/a
93n/a/* Random methods */
94n/a
95n/a
96n/a/* generates a random number on [0,0xffffffff]-interval */
97n/astatic uint32_t
98n/agenrand_int32(RandomObject *self)
99n/a{
100n/a uint32_t y;
101n/a static const uint32_t mag01[2] = {0x0U, MATRIX_A};
102n/a /* mag01[x] = x * MATRIX_A for x=0,1 */
103n/a uint32_t *mt;
104n/a
105n/a mt = self->state;
106n/a if (self->index >= N) { /* generate N words at one time */
107n/a int kk;
108n/a
109n/a for (kk=0;kk<N-M;kk++) {
110n/a y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
111n/a mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1U];
112n/a }
113n/a for (;kk<N-1;kk++) {
114n/a y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
115n/a mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1U];
116n/a }
117n/a y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
118n/a mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1U];
119n/a
120n/a self->index = 0;
121n/a }
122n/a
123n/a y = mt[self->index++];
124n/a y ^= (y >> 11);
125n/a y ^= (y << 7) & 0x9d2c5680U;
126n/a y ^= (y << 15) & 0xefc60000U;
127n/a y ^= (y >> 18);
128n/a return y;
129n/a}
130n/a
131n/a/* random_random is the function named genrand_res53 in the original code;
132n/a * generates a random number on [0,1) with 53-bit resolution; note that
133n/a * 9007199254740992 == 2**53; I assume they're spelling "/2**53" as
134n/a * multiply-by-reciprocal in the (likely vain) hope that the compiler will
135n/a * optimize the division away at compile-time. 67108864 is 2**26. In
136n/a * effect, a contains 27 random bits shifted left 26, and b fills in the
137n/a * lower 26 bits of the 53-bit numerator.
138n/a * The original code credited Isaku Wada for this algorithm, 2002/01/09.
139n/a */
140n/astatic PyObject *
141n/arandom_random(RandomObject *self)
142n/a{
143n/a uint32_t a=genrand_int32(self)>>5, b=genrand_int32(self)>>6;
144n/a return PyFloat_FromDouble((a*67108864.0+b)*(1.0/9007199254740992.0));
145n/a}
146n/a
147n/a/* initializes mt[N] with a seed */
148n/astatic void
149n/ainit_genrand(RandomObject *self, uint32_t s)
150n/a{
151n/a int mti;
152n/a uint32_t *mt;
153n/a
154n/a mt = self->state;
155n/a mt[0]= s;
156n/a for (mti=1; mti<N; mti++) {
157n/a mt[mti] =
158n/a (1812433253U * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
159n/a /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
160n/a /* In the previous versions, MSBs of the seed affect */
161n/a /* only MSBs of the array mt[]. */
162n/a /* 2002/01/09 modified by Makoto Matsumoto */
163n/a }
164n/a self->index = mti;
165n/a return;
166n/a}
167n/a
168n/a/* initialize by an array with array-length */
169n/a/* init_key is the array for initializing keys */
170n/a/* key_length is its length */
171n/astatic void
172n/ainit_by_array(RandomObject *self, uint32_t init_key[], size_t key_length)
173n/a{
174n/a size_t i, j, k; /* was signed in the original code. RDH 12/16/2002 */
175n/a uint32_t *mt;
176n/a
177n/a mt = self->state;
178n/a init_genrand(self, 19650218U);
179n/a i=1; j=0;
180n/a k = (N>key_length ? N : key_length);
181n/a for (; k; k--) {
182n/a mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525U))
183n/a + init_key[j] + (uint32_t)j; /* non linear */
184n/a i++; j++;
185n/a if (i>=N) { mt[0] = mt[N-1]; i=1; }
186n/a if (j>=key_length) j=0;
187n/a }
188n/a for (k=N-1; k; k--) {
189n/a mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941U))
190n/a - (uint32_t)i; /* non linear */
191n/a i++;
192n/a if (i>=N) { mt[0] = mt[N-1]; i=1; }
193n/a }
194n/a
195n/a mt[0] = 0x80000000U; /* MSB is 1; assuring non-zero initial array */
196n/a}
197n/a
198n/a/*
199n/a * The rest is Python-specific code, neither part of, nor derived from, the
200n/a * Twister download.
201n/a */
202n/a
203n/astatic int
204n/arandom_seed_urandom(RandomObject *self)
205n/a{
206n/a PY_UINT32_T key[N];
207n/a
208n/a if (_PyOS_URandomNonblock(key, sizeof(key)) < 0) {
209n/a return -1;
210n/a }
211n/a init_by_array(self, key, Py_ARRAY_LENGTH(key));
212n/a return 0;
213n/a}
214n/a
215n/astatic void
216n/arandom_seed_time_pid(RandomObject *self)
217n/a{
218n/a _PyTime_t now;
219n/a uint32_t key[5];
220n/a
221n/a now = _PyTime_GetSystemClock();
222n/a key[0] = (PY_UINT32_T)(now & 0xffffffffU);
223n/a key[1] = (PY_UINT32_T)(now >> 32);
224n/a
225n/a key[2] = (PY_UINT32_T)getpid();
226n/a
227n/a now = _PyTime_GetMonotonicClock();
228n/a key[3] = (PY_UINT32_T)(now & 0xffffffffU);
229n/a key[4] = (PY_UINT32_T)(now >> 32);
230n/a
231n/a init_by_array(self, key, Py_ARRAY_LENGTH(key));
232n/a}
233n/a
234n/astatic PyObject *
235n/arandom_seed(RandomObject *self, PyObject *args)
236n/a{
237n/a PyObject *result = NULL; /* guilty until proved innocent */
238n/a PyObject *n = NULL;
239n/a uint32_t *key = NULL;
240n/a size_t bits, keyused;
241n/a int res;
242n/a PyObject *arg = NULL;
243n/a
244n/a if (!PyArg_UnpackTuple(args, "seed", 0, 1, &arg))
245n/a return NULL;
246n/a
247n/a if (arg == NULL || arg == Py_None) {
248n/a if (random_seed_urandom(self) < 0) {
249n/a PyErr_Clear();
250n/a
251n/a /* Reading system entropy failed, fall back on the worst entropy:
252n/a use the current time and process identifier. */
253n/a random_seed_time_pid(self);
254n/a }
255n/a Py_RETURN_NONE;
256n/a }
257n/a
258n/a /* This algorithm relies on the number being unsigned.
259n/a * So: if the arg is a PyLong, use its absolute value.
260n/a * Otherwise use its hash value, cast to unsigned.
261n/a */
262n/a if (PyLong_Check(arg))
263n/a n = PyNumber_Absolute(arg);
264n/a else {
265n/a Py_hash_t hash = PyObject_Hash(arg);
266n/a if (hash == -1)
267n/a goto Done;
268n/a n = PyLong_FromSize_t((size_t)hash);
269n/a }
270n/a if (n == NULL)
271n/a goto Done;
272n/a
273n/a /* Now split n into 32-bit chunks, from the right. */
274n/a bits = _PyLong_NumBits(n);
275n/a if (bits == (size_t)-1 && PyErr_Occurred())
276n/a goto Done;
277n/a
278n/a /* Figure out how many 32-bit chunks this gives us. */
279n/a keyused = bits == 0 ? 1 : (bits - 1) / 32 + 1;
280n/a
281n/a /* Convert seed to byte sequence. */
282n/a key = (uint32_t *)PyMem_Malloc((size_t)4 * keyused);
283n/a if (key == NULL) {
284n/a PyErr_NoMemory();
285n/a goto Done;
286n/a }
287n/a res = _PyLong_AsByteArray((PyLongObject *)n,
288n/a (unsigned char *)key, keyused * 4,
289n/a PY_LITTLE_ENDIAN,
290n/a 0); /* unsigned */
291n/a if (res == -1) {
292n/a PyMem_Free(key);
293n/a goto Done;
294n/a }
295n/a
296n/a#if PY_BIG_ENDIAN
297n/a {
298n/a size_t i, j;
299n/a /* Reverse an array. */
300n/a for (i = 0, j = keyused - 1; i < j; i++, j--) {
301n/a uint32_t tmp = key[i];
302n/a key[i] = key[j];
303n/a key[j] = tmp;
304n/a }
305n/a }
306n/a#endif
307n/a init_by_array(self, key, keyused);
308n/a
309n/a Py_INCREF(Py_None);
310n/a result = Py_None;
311n/a
312n/aDone:
313n/a Py_XDECREF(n);
314n/a PyMem_Free(key);
315n/a return result;
316n/a}
317n/a
318n/astatic PyObject *
319n/arandom_getstate(RandomObject *self)
320n/a{
321n/a PyObject *state;
322n/a PyObject *element;
323n/a int i;
324n/a
325n/a state = PyTuple_New(N+1);
326n/a if (state == NULL)
327n/a return NULL;
328n/a for (i=0; i<N ; i++) {
329n/a element = PyLong_FromUnsignedLong(self->state[i]);
330n/a if (element == NULL)
331n/a goto Fail;
332n/a PyTuple_SET_ITEM(state, i, element);
333n/a }
334n/a element = PyLong_FromLong((long)(self->index));
335n/a if (element == NULL)
336n/a goto Fail;
337n/a PyTuple_SET_ITEM(state, i, element);
338n/a return state;
339n/a
340n/aFail:
341n/a Py_DECREF(state);
342n/a return NULL;
343n/a}
344n/a
345n/astatic PyObject *
346n/arandom_setstate(RandomObject *self, PyObject *state)
347n/a{
348n/a int i;
349n/a unsigned long element;
350n/a long index;
351n/a
352n/a if (!PyTuple_Check(state)) {
353n/a PyErr_SetString(PyExc_TypeError,
354n/a "state vector must be a tuple");
355n/a return NULL;
356n/a }
357n/a if (PyTuple_Size(state) != N+1) {
358n/a PyErr_SetString(PyExc_ValueError,
359n/a "state vector is the wrong size");
360n/a return NULL;
361n/a }
362n/a
363n/a for (i=0; i<N ; i++) {
364n/a element = PyLong_AsUnsignedLong(PyTuple_GET_ITEM(state, i));
365n/a if (element == (unsigned long)-1 && PyErr_Occurred())
366n/a return NULL;
367n/a self->state[i] = (uint32_t)element;
368n/a }
369n/a
370n/a index = PyLong_AsLong(PyTuple_GET_ITEM(state, i));
371n/a if (index == -1 && PyErr_Occurred())
372n/a return NULL;
373n/a if (index < 0 || index > N) {
374n/a PyErr_SetString(PyExc_ValueError, "invalid state");
375n/a return NULL;
376n/a }
377n/a self->index = (int)index;
378n/a
379n/a Py_RETURN_NONE;
380n/a}
381n/a
382n/astatic PyObject *
383n/arandom_getrandbits(RandomObject *self, PyObject *args)
384n/a{
385n/a int k, i, words;
386n/a uint32_t r;
387n/a uint32_t *wordarray;
388n/a PyObject *result;
389n/a
390n/a if (!PyArg_ParseTuple(args, "i:getrandbits", &k))
391n/a return NULL;
392n/a
393n/a if (k <= 0) {
394n/a PyErr_SetString(PyExc_ValueError,
395n/a "number of bits must be greater than zero");
396n/a return NULL;
397n/a }
398n/a
399n/a if (k <= 32) /* Fast path */
400n/a return PyLong_FromUnsignedLong(genrand_int32(self) >> (32 - k));
401n/a
402n/a words = (k - 1) / 32 + 1;
403n/a wordarray = (uint32_t *)PyMem_Malloc(words * 4);
404n/a if (wordarray == NULL) {
405n/a PyErr_NoMemory();
406n/a return NULL;
407n/a }
408n/a
409n/a /* Fill-out bits of long integer, by 32-bit words, from least significant
410n/a to most significant. */
411n/a#if PY_LITTLE_ENDIAN
412n/a for (i = 0; i < words; i++, k -= 32)
413n/a#else
414n/a for (i = words - 1; i >= 0; i--, k -= 32)
415n/a#endif
416n/a {
417n/a r = genrand_int32(self);
418n/a if (k < 32)
419n/a r >>= (32 - k); /* Drop least significant bits */
420n/a wordarray[i] = r;
421n/a }
422n/a
423n/a result = _PyLong_FromByteArray((unsigned char *)wordarray, words * 4,
424n/a PY_LITTLE_ENDIAN, 0 /* unsigned */);
425n/a PyMem_Free(wordarray);
426n/a return result;
427n/a}
428n/a
429n/astatic PyObject *
430n/arandom_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
431n/a{
432n/a RandomObject *self;
433n/a PyObject *tmp;
434n/a
435n/a if (type == &Random_Type && !_PyArg_NoKeywords("Random()", kwds))
436n/a return NULL;
437n/a
438n/a self = (RandomObject *)type->tp_alloc(type, 0);
439n/a if (self == NULL)
440n/a return NULL;
441n/a tmp = random_seed(self, args);
442n/a if (tmp == NULL) {
443n/a Py_DECREF(self);
444n/a return NULL;
445n/a }
446n/a Py_DECREF(tmp);
447n/a return (PyObject *)self;
448n/a}
449n/a
450n/astatic PyMethodDef random_methods[] = {
451n/a {"random", (PyCFunction)random_random, METH_NOARGS,
452n/a PyDoc_STR("random() -> x in the interval [0, 1).")},
453n/a {"seed", (PyCFunction)random_seed, METH_VARARGS,
454n/a PyDoc_STR("seed([n]) -> None. Defaults to current time.")},
455n/a {"getstate", (PyCFunction)random_getstate, METH_NOARGS,
456n/a PyDoc_STR("getstate() -> tuple containing the current state.")},
457n/a {"setstate", (PyCFunction)random_setstate, METH_O,
458n/a PyDoc_STR("setstate(state) -> None. Restores generator state.")},
459n/a {"getrandbits", (PyCFunction)random_getrandbits, METH_VARARGS,
460n/a PyDoc_STR("getrandbits(k) -> x. Generates an int with "
461n/a "k random bits.")},
462n/a {NULL, NULL} /* sentinel */
463n/a};
464n/a
465n/aPyDoc_STRVAR(random_doc,
466n/a"Random() -> create a random number generator with its own internal state.");
467n/a
468n/astatic PyTypeObject Random_Type = {
469n/a PyVarObject_HEAD_INIT(NULL, 0)
470n/a "_random.Random", /*tp_name*/
471n/a sizeof(RandomObject), /*tp_basicsize*/
472n/a 0, /*tp_itemsize*/
473n/a /* methods */
474n/a 0, /*tp_dealloc*/
475n/a 0, /*tp_print*/
476n/a 0, /*tp_getattr*/
477n/a 0, /*tp_setattr*/
478n/a 0, /*tp_reserved*/
479n/a 0, /*tp_repr*/
480n/a 0, /*tp_as_number*/
481n/a 0, /*tp_as_sequence*/
482n/a 0, /*tp_as_mapping*/
483n/a 0, /*tp_hash*/
484n/a 0, /*tp_call*/
485n/a 0, /*tp_str*/
486n/a PyObject_GenericGetAttr, /*tp_getattro*/
487n/a 0, /*tp_setattro*/
488n/a 0, /*tp_as_buffer*/
489n/a Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /*tp_flags*/
490n/a random_doc, /*tp_doc*/
491n/a 0, /*tp_traverse*/
492n/a 0, /*tp_clear*/
493n/a 0, /*tp_richcompare*/
494n/a 0, /*tp_weaklistoffset*/
495n/a 0, /*tp_iter*/
496n/a 0, /*tp_iternext*/
497n/a random_methods, /*tp_methods*/
498n/a 0, /*tp_members*/
499n/a 0, /*tp_getset*/
500n/a 0, /*tp_base*/
501n/a 0, /*tp_dict*/
502n/a 0, /*tp_descr_get*/
503n/a 0, /*tp_descr_set*/
504n/a 0, /*tp_dictoffset*/
505n/a 0, /*tp_init*/
506n/a 0, /*tp_alloc*/
507n/a random_new, /*tp_new*/
508n/a PyObject_Free, /*tp_free*/
509n/a 0, /*tp_is_gc*/
510n/a};
511n/a
512n/aPyDoc_STRVAR(module_doc,
513n/a"Module implements the Mersenne Twister random number generator.");
514n/a
515n/a
516n/astatic struct PyModuleDef _randommodule = {
517n/a PyModuleDef_HEAD_INIT,
518n/a "_random",
519n/a module_doc,
520n/a -1,
521n/a NULL,
522n/a NULL,
523n/a NULL,
524n/a NULL,
525n/a NULL
526n/a};
527n/a
528n/aPyMODINIT_FUNC
529n/aPyInit__random(void)
530n/a{
531n/a PyObject *m;
532n/a
533n/a if (PyType_Ready(&Random_Type) < 0)
534n/a return NULL;
535n/a m = PyModule_Create(&_randommodule);
536n/a if (m == NULL)
537n/a return NULL;
538n/a Py_INCREF(&Random_Type);
539n/a PyModule_AddObject(m, "Random", (PyObject *)&Random_Type);
540n/a return m;
541n/a}