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 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.keio.ac.jp/matumoto/emt.html 64 email: matumoto (at) math.keio.ac.jp 65 */ 66 67 /* ---------------------------------------------------------------*/ 68 69 #include "Python.h" 70 #include <time.h> /* for seeding to current time */ 71 72 /* Period parameters -- These are all magic. Don't change. */ 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 79 typedef struct { 80 PyObject_HEAD 81 unsigned long state[N]; 82 int index; 83 } RandomObject; 84 85 static PyTypeObject Random_Type; 86 87 #define RandomObject_Check(v) (Py_TYPE(v) == &Random_Type) 88 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_length) 172 { 173 unsigned int i, j, k; /* was signed in the original code. RDH 12/16/2002 */ 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 201 /* 202 * The rest is Python-specific code, neither part of, nor derived from, the 203 * Twister download. 204 */ 205 206 static PyObject * 207 random_seed(RandomObject *self, PyObject *args) 208 { 209 PyObject *result = NULL; /* guilty until proved innocent */ 210 PyObject *masklower = NULL; 211 PyObject *thirtytwo = NULL; 212 PyObject *n = NULL; 213 unsigned long *key = NULL; 214 unsigned long keymax; /* # of allocated slots in key */ 215 unsigned long keyused; /* # of used slots in key */ 216 int err; 217 218 PyObject *arg = NULL; 219 220 if (!PyArg_UnpackTuple(args, "seed", 0, 1, &arg)) 221 return NULL; 222 223 if (arg == NULL || arg == Py_None) { 224 time_t now; 225 226 time(&now); 227 init_genrand(self, (unsigned long)now); 228 Py_INCREF(Py_None); 229 return Py_None; 230 } 231 /* If the arg is an int or long, use its absolute value; else use 232 * the absolute value of its hash code. 233 */ 234 if (PyInt_Check(arg) || PyLong_Check(arg)) 235 n = PyNumber_Absolute(arg); 236 else { 237 long hash = PyObject_Hash(arg); 238 if (hash == -1) 239 goto Done; 240 n = PyLong_FromUnsignedLong((unsigned long)hash); 241 } 242 if (n == NULL) 243 goto Done; 244 245 /* Now split n into 32-bit chunks, from the right. Each piece is 246 * stored into key, which has a capacity of keymax chunks, of which 247 * keyused are filled. Alas, the repeated shifting makes this a 248 * quadratic-time algorithm; we'd really like to use 249 * _PyLong_AsByteArray here, but then we'd have to break into the 250 * long representation to figure out how big an array was needed 251 * in advance. 252 */ 253 keymax = 8; /* arbitrary; grows later if needed */ 254 keyused = 0; 255 key = (unsigned long *)PyMem_Malloc(keymax * sizeof(*key)); 256 if (key == NULL) 257 goto Done; 258 259 masklower = PyLong_FromUnsignedLong(0xffffffffU); 260 if (masklower == NULL) 261 goto Done; 262 thirtytwo = PyInt_FromLong(32L); 263 if (thirtytwo == NULL) 264 goto Done; 265 while ((err=PyObject_IsTrue(n))) { 266 PyObject *newn; 267 PyObject *pychunk; 268 unsigned long chunk; 269 270 if (err == -1) 271 goto Done; 272 pychunk = PyNumber_And(n, masklower); 273 if (pychunk == NULL) 274 goto Done; 275 chunk = PyLong_AsUnsignedLong(pychunk); 276 Py_DECREF(pychunk); 277 if (chunk == (unsigned long)-1 && PyErr_Occurred()) 278 goto Done; 279 newn = PyNumber_Rshift(n, thirtytwo); 280 if (newn == NULL) 281 goto Done; 282 Py_DECREF(n); 283 n = newn; 284 if (keyused >= keymax) { 285 unsigned long bigger = keymax << 1; 286 if ((bigger >> 1) != keymax) { 287 PyErr_NoMemory(); 288 goto Done; 289 } 290 key = (unsigned long *)PyMem_Realloc(key, 291 bigger * sizeof(*key)); 292 if (key == NULL) 293 goto Done; 294 keymax = bigger; 295 } 296 assert(keyused < keymax); 297 key[keyused++] = chunk; 298 } 299 300 if (keyused == 0) 301 key[keyused++] = 0UL; 302 result = init_by_array(self, key, keyused); 303 Done: 304 Py_XDECREF(masklower); 305 Py_XDECREF(thirtytwo); 306 Py_XDECREF(n); 307 PyMem_Free(key); 308 return result; 309 } 310 311 static PyObject * 312 random_getstate(RandomObject *self) 313 { 314 PyObject *state; 315 PyObject *element; 316 int i; 317 318 state = PyTuple_New(N+1); 319 if (state == NULL) 320 return NULL; 321 for (i=0; i<N ; i++) { 322 element = PyLong_FromUnsignedLong(self->state[i]); 323 if (element == NULL) 324 goto Fail; 325 PyTuple_SET_ITEM(state, i, element); 326 } 327 element = PyLong_FromLong((long)(self->index)); 328 if (element == NULL) 329 goto Fail; 330 PyTuple_SET_ITEM(state, i, element); 331 return state; 332 333 Fail: 334 Py_DECREF(state); 335 return NULL; 336 } 337 338 static PyObject * 339 random_setstate(RandomObject *self, PyObject *state) 340 { 341 int i; 342 unsigned long element; 343 long index; 344 345 if (!PyTuple_Check(state)) { 346 PyErr_SetString(PyExc_TypeError, 347 "state vector must be a tuple"); 348 return NULL; 349 } 350 if (PyTuple_Size(state) != N+1) { 351 PyErr_SetString(PyExc_ValueError, 352 "state vector is the wrong size"); 353 return NULL; 354 } 355 356 for (i=0; i<N ; i++) { 357 element = PyLong_AsUnsignedLong(PyTuple_GET_ITEM(state, i)); 358 if (element == (unsigned long)-1 && PyErr_Occurred()) 359 return NULL; 360 self->state[i] = element & 0xffffffffUL; /* Make sure we get sane state */ 361 } 362 363 index = PyLong_AsLong(PyTuple_GET_ITEM(state, i)); 364 if (index == -1 && PyErr_Occurred()) 365 return NULL; 366 self->index = (int)index; 367 368 Py_INCREF(Py_None); 369 return Py_None; 370 } 371 372 /* 373 Jumpahead should be a fast way advance the generator n-steps ahead, but 374 lacking a formula for that, the next best is to use n and the existing 375 state to create a new state far away from the original. 376 377 The generator uses constant spaced additive feedback, so shuffling the 378 state elements ought to produce a state which would not be encountered 379 (in the near term) by calls to random(). Shuffling is normally 380 implemented by swapping the ith element with another element ranging 381 from 0 to i inclusive. That allows the element to have the possibility 382 of not being moved. Since the goal is to produce a new, different 383 state, the swap element is ranged from 0 to i-1 inclusive. This assures 384 that each element gets moved at least once. 385 386 To make sure that consecutive calls to jumpahead(n) produce different 387 states (even in the rare case of involutory shuffles), i+1 is added to 388 each element at position i. Successive calls are then guaranteed to 389 have changing (growing) values as well as shuffled positions. 390 391 Finally, the self->index value is set to N so that the generator itself 392 kicks in on the next call to random(). This assures that all results 393 have been through the generator and do not just reflect alterations to 394 the underlying state. 395 */ 396 397 static PyObject * 398 random_jumpahead(RandomObject *self, PyObject *n) 399 { 400 long i, j; 401 PyObject *iobj; 402 PyObject *remobj; 403 unsigned long *mt, tmp, nonzero; 404 405 if (!PyInt_Check(n) && !PyLong_Check(n)) { 406 PyErr_Format(PyExc_TypeError, "jumpahead requires an " 407 "integer, not '%s'", 408 Py_TYPE(n)->tp_name); 409 return NULL; 410 } 411 412 mt = self->state; 413 for (i = N-1; i > 1; i--) { 414 iobj = PyInt_FromLong(i); 415 if (iobj == NULL) 416 return NULL; 417 remobj = PyNumber_Remainder(n, iobj); 418 Py_DECREF(iobj); 419 if (remobj == NULL) 420 return NULL; 421 j = PyInt_AsLong(remobj); 422 Py_DECREF(remobj); 423 if (j == -1L && PyErr_Occurred()) 424 return NULL; 425 tmp = mt[i]; 426 mt[i] = mt[j]; 427 mt[j] = tmp; 428 } 429 430 nonzero = 0; 431 for (i = 1; i < N; i++) { 432 mt[i] += i+1; 433 mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ 434 nonzero |= mt[i]; 435 } 436 437 /* Ensure the state is nonzero: in the unlikely event that mt[1] through 438 mt[N-1] are all zero, set the MSB of mt[0] (see issue #14591). In the 439 normal case, we fall back to the pre-issue 14591 behaviour for mt[0]. */ 440 if (nonzero) { 441 mt[0] += 1; 442 mt[0] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ 443 } 444 else { 445 mt[0] = 0x80000000UL; 446 } 447 448 self->index = N; 449 Py_INCREF(Py_None); 450 return Py_None; 451 } 452 453 static PyObject * 454 random_getrandbits(RandomObject *self, PyObject *args) 455 { 456 int k, i, bytes; 457 unsigned long r; 458 unsigned char *bytearray; 459 PyObject *result; 460 461 if (!PyArg_ParseTuple(args, "i:getrandbits", &k)) 462 return NULL; 463 464 if (k <= 0) { 465 PyErr_SetString(PyExc_ValueError, 466 "number of bits must be greater than zero"); 467 return NULL; 468 } 469 470 bytes = ((k - 1) / 32 + 1) * 4; 471 bytearray = (unsigned char *)PyMem_Malloc(bytes); 472 if (bytearray == NULL) { 473 PyErr_NoMemory(); 474 return NULL; 475 } 476 477 /* Fill-out whole words, byte-by-byte to avoid endianness issues */ 478 for (i=0 ; i<bytes ; i+=4, k-=32) { 479 r = genrand_int32(self); 480 if (k < 32) 481 r >>= (32 - k); 482 bytearray[i+0] = (unsigned char)r; 483 bytearray[i+1] = (unsigned char)(r >> 8); 484 bytearray[i+2] = (unsigned char)(r >> 16); 485 bytearray[i+3] = (unsigned char)(r >> 24); 486 } 487 488 /* little endian order to match bytearray assignment order */ 489 result = _PyLong_FromByteArray(bytearray, bytes, 1, 0); 490 PyMem_Free(bytearray); 491 return result; 492 } 493 494 static PyObject * 495 random_new(PyTypeObject *type, PyObject *args, PyObject *kwds) 496 { 497 RandomObject *self; 498 PyObject *tmp; 499 500 if (type == &Random_Type && !_PyArg_NoKeywords("Random()", kwds)) 501 return NULL; 502 503 self = (RandomObject *)type->tp_alloc(type, 0); 504 if (self == NULL) 505 return NULL; 506 tmp = random_seed(self, args); 507 if (tmp == NULL) { 508 Py_DECREF(self); 509 return NULL; 510 } 511 Py_DECREF(tmp); 512 return (PyObject *)self; 513 } 514 515 static PyMethodDef random_methods[] = { 516 {"random", (PyCFunction)random_random, METH_NOARGS, 517 PyDoc_STR("random() -> x in the interval [0, 1).")}, 518 {"seed", (PyCFunction)random_seed, METH_VARARGS, 519 PyDoc_STR("seed([n]) -> None. Defaults to current time.")}, 520 {"getstate", (PyCFunction)random_getstate, METH_NOARGS, 521 PyDoc_STR("getstate() -> tuple containing the current state.")}, 522 {"setstate", (PyCFunction)random_setstate, METH_O, 523 PyDoc_STR("setstate(state) -> None. Restores generator state.")}, 524 {"jumpahead", (PyCFunction)random_jumpahead, METH_O, 525 PyDoc_STR("jumpahead(int) -> None. Create new state from " 526 "existing state and integer.")}, 527 {"getrandbits", (PyCFunction)random_getrandbits, METH_VARARGS, 528 PyDoc_STR("getrandbits(k) -> x. Generates a long int with " 529 "k random bits.")}, 530 {NULL, NULL} /* sentinel */ 531 }; 532 533 PyDoc_STRVAR(random_doc, 534 "Random() -> create a random number generator with its own internal state."); 535 536 static PyTypeObject Random_Type = { 537 PyVarObject_HEAD_INIT(NULL, 0) 538 "_random.Random", /*tp_name*/ 539 sizeof(RandomObject), /*tp_basicsize*/ 540 0, /*tp_itemsize*/ 541 /* methods */ 542 0, /*tp_dealloc*/ 543 0, /*tp_print*/ 544 0, /*tp_getattr*/ 545 0, /*tp_setattr*/ 546 0, /*tp_compare*/ 547 0, /*tp_repr*/ 548 0, /*tp_as_number*/ 549 0, /*tp_as_sequence*/ 550 0, /*tp_as_mapping*/ 551 0, /*tp_hash*/ 552 0, /*tp_call*/ 553 0, /*tp_str*/ 554 PyObject_GenericGetAttr, /*tp_getattro*/ 555 0, /*tp_setattro*/ 556 0, /*tp_as_buffer*/ 557 Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /*tp_flags*/ 558 random_doc, /*tp_doc*/ 559 0, /*tp_traverse*/ 560 0, /*tp_clear*/ 561 0, /*tp_richcompare*/ 562 0, /*tp_weaklistoffset*/ 563 0, /*tp_iter*/ 564 0, /*tp_iternext*/ 565 random_methods, /*tp_methods*/ 566 0, /*tp_members*/ 567 0, /*tp_getset*/ 568 0, /*tp_base*/ 569 0, /*tp_dict*/ 570 0, /*tp_descr_get*/ 571 0, /*tp_descr_set*/ 572 0, /*tp_dictoffset*/ 573 0, /*tp_init*/ 574 0, /*tp_alloc*/ 575 random_new, /*tp_new*/ 576 _PyObject_Del, /*tp_free*/ 577 0, /*tp_is_gc*/ 578 }; 579 580 PyDoc_STRVAR(module_doc, 581 "Module implements the Mersenne Twister random number generator."); 582 583 PyMODINIT_FUNC 584 init_random(void) 585 { 586 PyObject *m; 587 588 if (PyType_Ready(&Random_Type) < 0) 589 return; 590 m = Py_InitModule3("_random", NULL, module_doc); 591 if (m == NULL) 592 return; 593 Py_INCREF(&Random_Type); 594 PyModule_AddObject(m, "Random", (PyObject *)&Random_Type); 595 } 596
