Rietveld Code Review Tool
Help | Bug tracker | Discussion group | Source code | Sign in
(411)

Side by Side Diff: Modules/_randommodule.c

Issue 13704: Random number generator in Python core
Patch Set: Created 8 years ago
Left:
Right:
Use n/p to move between diff chunks; N/P to move between comments. Please Sign in to add in-line comments.
Jump to:
View unified diff | Download patch
OLDNEW
1 /* Random objects */ 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 O R
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" 2 #include "Python.h"
3 #include "pyrandom.h"
70 #include <time.h> /* for seeding to current time */ 4 #include <time.h> /* for seeding to current time */
71 5
72 /* Period parameters -- These are all magic. Don't change. */ 6 #define N _Py_MT_N
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 7
79 typedef struct { 8 typedef struct {
80 PyObject_HEAD 9 PyObject_HEAD
81 unsigned long state[N]; 10 _Py_MT_RandomState rs;
82 int index;
83 } RandomObject; 11 } RandomObject;
84 12
85 static PyTypeObject Random_Type; 13 static PyTypeObject Random_Type;
86 14
87 #define RandomObject_Check(v) (Py_TYPE(v) == &Random_Type) 15 #define RandomObject_Check(v) (Py_TYPE(v) == &Random_Type)
88 16
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_le ngth)
172 {
173 unsigned int i, j, k; /* was signed in the original code. RDH 12/16/20 02 */
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 17
201 /* 18 /*
202 * The rest is Python-specific code, neither part of, nor derived from, the 19 * The rest is Python-specific code, neither part of, nor derived from, the
203 * Twister download. 20 * Twister download.
204 */ 21 */
205 22
206 static PyObject * 23 static PyObject *
24 random_random(RandomObject *self)
25 {
26 return PyFloat_FromDouble(_Py_MT_GenRand_Res53(&self->rs));
27 }
28
29
30 static PyObject *
207 random_seed(RandomObject *self, PyObject *args) 31 random_seed(RandomObject *self, PyObject *args)
208 { 32 {
209 PyObject *result = NULL; /* guilty until proved innocent */ 33 PyObject *result = NULL; /* guilty until proved innocent */
210 PyObject *masklower = NULL; 34 PyObject *masklower = NULL;
211 PyObject *thirtytwo = NULL; 35 PyObject *thirtytwo = NULL;
212 PyObject *n = NULL; 36 PyObject *n = NULL;
213 unsigned long *key = NULL; 37 unsigned long *key = NULL;
214 unsigned long keymax; /* # of allocated slots in key */ 38 unsigned long keymax; /* # of allocated slots in key */
215 unsigned long keyused; /* # of used slots in key */ 39 unsigned long keyused; /* # of used slots in key */
216 int err; 40 int err;
217 41
218 PyObject *arg = NULL; 42 PyObject *arg = NULL;
219 43
220 if (!PyArg_UnpackTuple(args, "seed", 0, 1, &arg)) 44 if (!PyArg_UnpackTuple(args, "seed", 0, 1, &arg))
221 return NULL; 45 return NULL;
222 46
223 if (arg == NULL || arg == Py_None) { 47 if (arg == NULL || arg == Py_None) {
224 time_t now; 48 time_t now;
225 49
226 time(&now); 50 time(&now);
227 init_genrand(self, (unsigned long)now); 51 _Py_MT_GenRand_Init(&self->rs, (unsigned long)now);
228 Py_INCREF(Py_None); 52 Py_INCREF(Py_None);
229 return Py_None; 53 return Py_None;
230 } 54 }
231 /* If the arg is an int or long, use its absolute value; else use 55 /* If the arg is an int or long, use its absolute value; else use
232 * the absolute value of its hash code. 56 * the absolute value of its hash code.
233 */ 57 */
234 if (PyLong_Check(arg)) 58 if (PyLong_Check(arg))
235 n = PyNumber_Absolute(arg); 59 n = PyNumber_Absolute(arg);
236 else { 60 else {
237 long hash = PyObject_Hash(arg); 61 long hash = PyObject_Hash(arg);
(...skipping 54 matching lines...) Expand 10 before | Expand all | Expand 10 after
292 if (key == NULL) 116 if (key == NULL)
293 goto Done; 117 goto Done;
294 keymax = bigger; 118 keymax = bigger;
295 } 119 }
296 assert(keyused < keymax); 120 assert(keyused < keymax);
297 key[keyused++] = chunk; 121 key[keyused++] = chunk;
298 } 122 }
299 123
300 if (keyused == 0) 124 if (keyused == 0)
301 key[keyused++] = 0UL; 125 key[keyused++] = 0UL;
302 result = init_by_array(self, key, keyused); 126 _Py_MT_GenRand_InitArray(&self->rs, key, keyused);
127 Py_INCREF(Py_None);
128 result = Py_None;
303 Done: 129 Done:
304 Py_XDECREF(masklower); 130 Py_XDECREF(masklower);
305 Py_XDECREF(thirtytwo); 131 Py_XDECREF(thirtytwo);
306 Py_XDECREF(n); 132 Py_XDECREF(n);
307 PyMem_Free(key); 133 PyMem_Free(key);
308 return result; 134 return result;
309 } 135 }
310 136
311 static PyObject * 137 static PyObject *
312 random_getstate(RandomObject *self) 138 random_getstate(RandomObject *self)
313 { 139 {
314 PyObject *state; 140 PyObject *state;
315 PyObject *element; 141 PyObject *element;
316 int i; 142 int i;
317 143
318 state = PyTuple_New(N+1); 144 state = PyTuple_New(N+1);
319 if (state == NULL) 145 if (state == NULL)
320 return NULL; 146 return NULL;
321 for (i=0; i<N ; i++) { 147 for (i=0; i<N ; i++) {
322 element = PyLong_FromUnsignedLong(self->state[i]); 148 element = PyLong_FromUnsignedLong(self->rs.state[i]);
323 if (element == NULL) 149 if (element == NULL)
324 goto Fail; 150 goto Fail;
325 PyTuple_SET_ITEM(state, i, element); 151 PyTuple_SET_ITEM(state, i, element);
326 } 152 }
327 element = PyLong_FromLong((long)(self->index)); 153 element = PyLong_FromLong((long)(self->rs.index));
328 if (element == NULL) 154 if (element == NULL)
329 goto Fail; 155 goto Fail;
330 PyTuple_SET_ITEM(state, i, element); 156 PyTuple_SET_ITEM(state, i, element);
331 return state; 157 return state;
332 158
333 Fail: 159 Fail:
334 Py_DECREF(state); 160 Py_DECREF(state);
335 return NULL; 161 return NULL;
336 } 162 }
337 163
(...skipping 12 matching lines...) Expand all
350 if (PyTuple_Size(state) != N+1) { 176 if (PyTuple_Size(state) != N+1) {
351 PyErr_SetString(PyExc_ValueError, 177 PyErr_SetString(PyExc_ValueError,
352 "state vector is the wrong size"); 178 "state vector is the wrong size");
353 return NULL; 179 return NULL;
354 } 180 }
355 181
356 for (i=0; i<N ; i++) { 182 for (i=0; i<N ; i++) {
357 element = PyLong_AsUnsignedLong(PyTuple_GET_ITEM(state, i)); 183 element = PyLong_AsUnsignedLong(PyTuple_GET_ITEM(state, i));
358 if (element == (unsigned long)-1 && PyErr_Occurred()) 184 if (element == (unsigned long)-1 && PyErr_Occurred())
359 return NULL; 185 return NULL;
360 self->state[i] = element & 0xffffffffUL; /* Make sure we get sane state */ 186 self->rs.state[i] = element & 0xffffffffUL; /* Make sure we get sane sta te */
361 } 187 }
362 188
363 index = PyLong_AsLong(PyTuple_GET_ITEM(state, i)); 189 index = PyLong_AsLong(PyTuple_GET_ITEM(state, i));
364 if (index == -1 && PyErr_Occurred()) 190 if (index == -1 && PyErr_Occurred())
365 return NULL; 191 return NULL;
366 self->index = (int)index; 192 self->rs.index = (int)index;
367 193
368 Py_INCREF(Py_None); 194 Py_INCREF(Py_None);
369 return Py_None; 195 return Py_None;
370 } 196 }
371 197
372 static PyObject * 198 static PyObject *
373 random_getrandbits(RandomObject *self, PyObject *args) 199 random_getrandbits(RandomObject *self, PyObject *args)
374 { 200 {
375 int k, i, bytes; 201 int k, i, bytes;
376 unsigned long r; 202 unsigned long r;
(...skipping 11 matching lines...) Expand all
388 214
389 bytes = ((k - 1) / 32 + 1) * 4; 215 bytes = ((k - 1) / 32 + 1) * 4;
390 bytearray = (unsigned char *)PyMem_Malloc(bytes); 216 bytearray = (unsigned char *)PyMem_Malloc(bytes);
391 if (bytearray == NULL) { 217 if (bytearray == NULL) {
392 PyErr_NoMemory(); 218 PyErr_NoMemory();
393 return NULL; 219 return NULL;
394 } 220 }
395 221
396 /* Fill-out whole words, byte-by-byte to avoid endianness issues */ 222 /* Fill-out whole words, byte-by-byte to avoid endianness issues */
397 for (i=0 ; i<bytes ; i+=4, k-=32) { 223 for (i=0 ; i<bytes ; i+=4, k-=32) {
398 r = genrand_int32(self); 224 r = _Py_MT_GenRand_Int32(&self->rs);
399 if (k < 32) 225 if (k < 32)
400 r >>= (32 - k); 226 r >>= (32 - k);
401 bytearray[i+0] = (unsigned char)r; 227 bytearray[i+0] = (unsigned char)r;
402 bytearray[i+1] = (unsigned char)(r >> 8); 228 bytearray[i+1] = (unsigned char)(r >> 8);
403 bytearray[i+2] = (unsigned char)(r >> 16); 229 bytearray[i+2] = (unsigned char)(r >> 16);
404 bytearray[i+3] = (unsigned char)(r >> 24); 230 bytearray[i+3] = (unsigned char)(r >> 24);
405 } 231 }
406 232
407 /* little endian order to match bytearray assignment order */ 233 /* little endian order to match bytearray assignment order */
408 result = _PyLong_FromByteArray(bytearray, bytes, 1, 0); 234 result = _PyLong_FromByteArray(bytearray, bytes, 1, 0);
(...skipping 107 matching lines...) Expand 10 before | Expand all | Expand 10 after
516 342
517 if (PyType_Ready(&Random_Type) < 0) 343 if (PyType_Ready(&Random_Type) < 0)
518 return NULL; 344 return NULL;
519 m = PyModule_Create(&_randommodule); 345 m = PyModule_Create(&_randommodule);
520 if (m == NULL) 346 if (m == NULL)
521 return NULL; 347 return NULL;
522 Py_INCREF(&Random_Type); 348 Py_INCREF(&Random_Type);
523 PyModule_AddObject(m, "Random", (PyObject *)&Random_Type); 349 PyModule_AddObject(m, "Random", (PyObject *)&Random_Type);
524 return m; 350 return m;
525 } 351 }
OLDNEW

RSS Feeds Recent Issues | This issue
This is Rietveld 894c83f36cb7+