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@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 79typedef struct { 80 PyObject_HEAD 81 unsigned long state[N]; 82 int index; 83} RandomObject; 84 85static 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 */ 94static unsigned long 95genrand_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 */ 137static PyObject * 138random_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 */ 145static void 146init_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 */ 170static PyObject * 171init_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 206static PyObject * 207random_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); 303Done: 304 Py_XDECREF(masklower); 305 Py_XDECREF(thirtytwo); 306 Py_XDECREF(n); 307 PyMem_Free(key); 308 return result; 309} 310 311static PyObject * 312random_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 333Fail: 334 Py_DECREF(state); 335 return NULL; 336} 337 338static PyObject * 339random_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 if (index < 0 || index > N) { 367 PyErr_SetString(PyExc_ValueError, "invalid state"); 368 return NULL; 369 } 370 self->index = (int)index; 371 372 Py_INCREF(Py_None); 373 return Py_None; 374} 375 376/* 377Jumpahead should be a fast way advance the generator n-steps ahead, but 378lacking a formula for that, the next best is to use n and the existing 379state to create a new state far away from the original. 380 381The generator uses constant spaced additive feedback, so shuffling the 382state elements ought to produce a state which would not be encountered 383(in the near term) by calls to random(). Shuffling is normally 384implemented by swapping the ith element with another element ranging 385from 0 to i inclusive. That allows the element to have the possibility 386of not being moved. Since the goal is to produce a new, different 387state, the swap element is ranged from 0 to i-1 inclusive. This assures 388that each element gets moved at least once. 389 390To make sure that consecutive calls to jumpahead(n) produce different 391states (even in the rare case of involutory shuffles), i+1 is added to 392each element at position i. Successive calls are then guaranteed to 393have changing (growing) values as well as shuffled positions. 394 395Finally, the self->index value is set to N so that the generator itself 396kicks in on the next call to random(). This assures that all results 397have been through the generator and do not just reflect alterations to 398the underlying state. 399*/ 400 401static PyObject * 402random_jumpahead(RandomObject *self, PyObject *n) 403{ 404 long i, j; 405 PyObject *iobj; 406 PyObject *remobj; 407 unsigned long *mt, tmp, nonzero; 408 409 if (!PyInt_Check(n) && !PyLong_Check(n)) { 410 PyErr_Format(PyExc_TypeError, "jumpahead requires an " 411 "integer, not '%s'", 412 Py_TYPE(n)->tp_name); 413 return NULL; 414 } 415 416 mt = self->state; 417 for (i = N-1; i > 1; i--) { 418 iobj = PyInt_FromLong(i); 419 if (iobj == NULL) 420 return NULL; 421 remobj = PyNumber_Remainder(n, iobj); 422 Py_DECREF(iobj); 423 if (remobj == NULL) 424 return NULL; 425 j = PyInt_AsLong(remobj); 426 Py_DECREF(remobj); 427 if (j == -1L && PyErr_Occurred()) 428 return NULL; 429 tmp = mt[i]; 430 mt[i] = mt[j]; 431 mt[j] = tmp; 432 } 433 434 nonzero = 0; 435 for (i = 1; i < N; i++) { 436 mt[i] += i+1; 437 mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ 438 nonzero |= mt[i]; 439 } 440 441 /* Ensure the state is nonzero: in the unlikely event that mt[1] through 442 mt[N-1] are all zero, set the MSB of mt[0] (see issue #14591). In the 443 normal case, we fall back to the pre-issue 14591 behaviour for mt[0]. */ 444 if (nonzero) { 445 mt[0] += 1; 446 mt[0] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ 447 } 448 else { 449 mt[0] = 0x80000000UL; 450 } 451 452 self->index = N; 453 Py_INCREF(Py_None); 454 return Py_None; 455} 456 457static PyObject * 458random_getrandbits(RandomObject *self, PyObject *args) 459{ 460 int k, i, bytes; 461 unsigned long r; 462 unsigned char *bytearray; 463 PyObject *result; 464 465 if (!PyArg_ParseTuple(args, "i:getrandbits", &k)) 466 return NULL; 467 468 if (k <= 0) { 469 PyErr_SetString(PyExc_ValueError, 470 "number of bits must be greater than zero"); 471 return NULL; 472 } 473 474 bytes = ((k - 1) / 32 + 1) * 4; 475 bytearray = (unsigned char *)PyMem_Malloc(bytes); 476 if (bytearray == NULL) { 477 PyErr_NoMemory(); 478 return NULL; 479 } 480 481 /* Fill-out whole words, byte-by-byte to avoid endianness issues */ 482 for (i=0 ; i<bytes ; i+=4, k-=32) { 483 r = genrand_int32(self); 484 if (k < 32) 485 r >>= (32 - k); 486 bytearray[i+0] = (unsigned char)r; 487 bytearray[i+1] = (unsigned char)(r >> 8); 488 bytearray[i+2] = (unsigned char)(r >> 16); 489 bytearray[i+3] = (unsigned char)(r >> 24); 490 } 491 492 /* little endian order to match bytearray assignment order */ 493 result = _PyLong_FromByteArray(bytearray, bytes, 1, 0); 494 PyMem_Free(bytearray); 495 return result; 496} 497 498static PyObject * 499random_new(PyTypeObject *type, PyObject *args, PyObject *kwds) 500{ 501 RandomObject *self; 502 PyObject *tmp; 503 504 if (type == &Random_Type && !_PyArg_NoKeywords("Random()", kwds)) 505 return NULL; 506 507 self = (RandomObject *)type->tp_alloc(type, 0); 508 if (self == NULL) 509 return NULL; 510 tmp = random_seed(self, args); 511 if (tmp == NULL) { 512 Py_DECREF(self); 513 return NULL; 514 } 515 Py_DECREF(tmp); 516 return (PyObject *)self; 517} 518 519static PyMethodDef random_methods[] = { 520 {"random", (PyCFunction)random_random, METH_NOARGS, 521 PyDoc_STR("random() -> x in the interval [0, 1).")}, 522 {"seed", (PyCFunction)random_seed, METH_VARARGS, 523 PyDoc_STR("seed([n]) -> None. Defaults to current time.")}, 524 {"getstate", (PyCFunction)random_getstate, METH_NOARGS, 525 PyDoc_STR("getstate() -> tuple containing the current state.")}, 526 {"setstate", (PyCFunction)random_setstate, METH_O, 527 PyDoc_STR("setstate(state) -> None. Restores generator state.")}, 528 {"jumpahead", (PyCFunction)random_jumpahead, METH_O, 529 PyDoc_STR("jumpahead(int) -> None. Create new state from " 530 "existing state and integer.")}, 531 {"getrandbits", (PyCFunction)random_getrandbits, METH_VARARGS, 532 PyDoc_STR("getrandbits(k) -> x. Generates a long int with " 533 "k random bits.")}, 534 {NULL, NULL} /* sentinel */ 535}; 536 537PyDoc_STRVAR(random_doc, 538"Random() -> create a random number generator with its own internal state."); 539 540static PyTypeObject Random_Type = { 541 PyVarObject_HEAD_INIT(NULL, 0) 542 "_random.Random", /*tp_name*/ 543 sizeof(RandomObject), /*tp_basicsize*/ 544 0, /*tp_itemsize*/ 545 /* methods */ 546 0, /*tp_dealloc*/ 547 0, /*tp_print*/ 548 0, /*tp_getattr*/ 549 0, /*tp_setattr*/ 550 0, /*tp_compare*/ 551 0, /*tp_repr*/ 552 0, /*tp_as_number*/ 553 0, /*tp_as_sequence*/ 554 0, /*tp_as_mapping*/ 555 0, /*tp_hash*/ 556 0, /*tp_call*/ 557 0, /*tp_str*/ 558 PyObject_GenericGetAttr, /*tp_getattro*/ 559 0, /*tp_setattro*/ 560 0, /*tp_as_buffer*/ 561 Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /*tp_flags*/ 562 random_doc, /*tp_doc*/ 563 0, /*tp_traverse*/ 564 0, /*tp_clear*/ 565 0, /*tp_richcompare*/ 566 0, /*tp_weaklistoffset*/ 567 0, /*tp_iter*/ 568 0, /*tp_iternext*/ 569 random_methods, /*tp_methods*/ 570 0, /*tp_members*/ 571 0, /*tp_getset*/ 572 0, /*tp_base*/ 573 0, /*tp_dict*/ 574 0, /*tp_descr_get*/ 575 0, /*tp_descr_set*/ 576 0, /*tp_dictoffset*/ 577 0, /*tp_init*/ 578 0, /*tp_alloc*/ 579 random_new, /*tp_new*/ 580 _PyObject_Del, /*tp_free*/ 581 0, /*tp_is_gc*/ 582}; 583 584PyDoc_STRVAR(module_doc, 585"Module implements the Mersenne Twister random number generator."); 586 587PyMODINIT_FUNC 588init_random(void) 589{ 590 PyObject *m; 591 592 if (PyType_Ready(&Random_Type) < 0) 593 return; 594 m = Py_InitModule3("_random", NULL, module_doc); 595 if (m == NULL) 596 return; 597 Py_INCREF(&Random_Type); 598 PyModule_AddObject(m, "Random", (PyObject *)&Random_Type); 599} 600