/Modules/_randommodule.c
C | 580 lines | 391 code | 60 blank | 129 comment | 89 complexity | f79f21b2c33a3f4dda811a0cc29d6f17 MD5 | raw file
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 except for tabbing. 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 self->index = (int)index; 367 368 Py_INCREF(Py_None); 369 return Py_None; 370} 371 372/* 373Jumpahead should be a fast way advance the generator n-steps ahead, but 374lacking a formula for that, the next best is to use n and the existing 375state to create a new state far away from the original. 376 377The generator uses constant spaced additive feedback, so shuffling the 378state elements ought to produce a state which would not be encountered 379(in the near term) by calls to random(). Shuffling is normally 380implemented by swapping the ith element with another element ranging 381from 0 to i inclusive. That allows the element to have the possibility 382of not being moved. Since the goal is to produce a new, different 383state, the swap element is ranged from 0 to i-1 inclusive. This assures 384that each element gets moved at least once. 385 386To make sure that consecutive calls to jumpahead(n) produce different 387states (even in the rare case of involutory shuffles), i+1 is added to 388each element at position i. Successive calls are then guaranteed to 389have changing (growing) values as well as shuffled positions. 390 391Finally, the self->index value is set to N so that the generator itself 392kicks in on the next call to random(). This assures that all results 393have been through the generator and do not just reflect alterations to 394the underlying state. 395*/ 396 397static PyObject * 398random_jumpahead(RandomObject *self, PyObject *n) 399{ 400 long i, j; 401 PyObject *iobj; 402 PyObject *remobj; 403 unsigned long *mt, tmp; 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 for (i = 0; i < N; i++) 431 mt[i] += i+1; 432 433 self->index = N; 434 Py_INCREF(Py_None); 435 return Py_None; 436} 437 438static PyObject * 439random_getrandbits(RandomObject *self, PyObject *args) 440{ 441 int k, i, bytes; 442 unsigned long r; 443 unsigned char *bytearray; 444 PyObject *result; 445 446 if (!PyArg_ParseTuple(args, "i:getrandbits", &k)) 447 return NULL; 448 449 if (k <= 0) { 450 PyErr_SetString(PyExc_ValueError, 451 "number of bits must be greater than zero"); 452 return NULL; 453 } 454 455 bytes = ((k - 1) / 32 + 1) * 4; 456 bytearray = (unsigned char *)PyMem_Malloc(bytes); 457 if (bytearray == NULL) { 458 PyErr_NoMemory(); 459 return NULL; 460 } 461 462 /* Fill-out whole words, byte-by-byte to avoid endianness issues */ 463 for (i=0 ; i<bytes ; i+=4, k-=32) { 464 r = genrand_int32(self); 465 if (k < 32) 466 r >>= (32 - k); 467 bytearray[i+0] = (unsigned char)r; 468 bytearray[i+1] = (unsigned char)(r >> 8); 469 bytearray[i+2] = (unsigned char)(r >> 16); 470 bytearray[i+3] = (unsigned char)(r >> 24); 471 } 472 473 /* little endian order to match bytearray assignment order */ 474 result = _PyLong_FromByteArray(bytearray, bytes, 1, 0); 475 PyMem_Free(bytearray); 476 return result; 477} 478 479static PyObject * 480random_new(PyTypeObject *type, PyObject *args, PyObject *kwds) 481{ 482 RandomObject *self; 483 PyObject *tmp; 484 485 if (type == &Random_Type && !_PyArg_NoKeywords("Random()", kwds)) 486 return NULL; 487 488 self = (RandomObject *)type->tp_alloc(type, 0); 489 if (self == NULL) 490 return NULL; 491 tmp = random_seed(self, args); 492 if (tmp == NULL) { 493 Py_DECREF(self); 494 return NULL; 495 } 496 Py_DECREF(tmp); 497 return (PyObject *)self; 498} 499 500static PyMethodDef random_methods[] = { 501 {"random", (PyCFunction)random_random, METH_NOARGS, 502 PyDoc_STR("random() -> x in the interval [0, 1).")}, 503 {"seed", (PyCFunction)random_seed, METH_VARARGS, 504 PyDoc_STR("seed([n]) -> None. Defaults to current time.")}, 505 {"getstate", (PyCFunction)random_getstate, METH_NOARGS, 506 PyDoc_STR("getstate() -> tuple containing the current state.")}, 507 {"setstate", (PyCFunction)random_setstate, METH_O, 508 PyDoc_STR("setstate(state) -> None. Restores generator state.")}, 509 {"jumpahead", (PyCFunction)random_jumpahead, METH_O, 510 PyDoc_STR("jumpahead(int) -> None. Create new state from " 511 "existing state and integer.")}, 512 {"getrandbits", (PyCFunction)random_getrandbits, METH_VARARGS, 513 PyDoc_STR("getrandbits(k) -> x. Generates a long int with " 514 "k random bits.")}, 515 {NULL, NULL} /* sentinel */ 516}; 517 518PyDoc_STRVAR(random_doc, 519"Random() -> create a random number generator with its own internal state."); 520 521static PyTypeObject Random_Type = { 522 PyVarObject_HEAD_INIT(NULL, 0) 523 "_random.Random", /*tp_name*/ 524 sizeof(RandomObject), /*tp_basicsize*/ 525 0, /*tp_itemsize*/ 526 /* methods */ 527 0, /*tp_dealloc*/ 528 0, /*tp_print*/ 529 0, /*tp_getattr*/ 530 0, /*tp_setattr*/ 531 0, /*tp_compare*/ 532 0, /*tp_repr*/ 533 0, /*tp_as_number*/ 534 0, /*tp_as_sequence*/ 535 0, /*tp_as_mapping*/ 536 0, /*tp_hash*/ 537 0, /*tp_call*/ 538 0, /*tp_str*/ 539 PyObject_GenericGetAttr, /*tp_getattro*/ 540 0, /*tp_setattro*/ 541 0, /*tp_as_buffer*/ 542 Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /*tp_flags*/ 543 random_doc, /*tp_doc*/ 544 0, /*tp_traverse*/ 545 0, /*tp_clear*/ 546 0, /*tp_richcompare*/ 547 0, /*tp_weaklistoffset*/ 548 0, /*tp_iter*/ 549 0, /*tp_iternext*/ 550 random_methods, /*tp_methods*/ 551 0, /*tp_members*/ 552 0, /*tp_getset*/ 553 0, /*tp_base*/ 554 0, /*tp_dict*/ 555 0, /*tp_descr_get*/ 556 0, /*tp_descr_set*/ 557 0, /*tp_dictoffset*/ 558 0, /*tp_init*/ 559 0, /*tp_alloc*/ 560 random_new, /*tp_new*/ 561 _PyObject_Del, /*tp_free*/ 562 0, /*tp_is_gc*/ 563}; 564 565PyDoc_STRVAR(module_doc, 566"Module implements the Mersenne Twister random number generator."); 567 568PyMODINIT_FUNC 569init_random(void) 570{ 571 PyObject *m; 572 573 if (PyType_Ready(&Random_Type) < 0) 574 return; 575 m = Py_InitModule3("_random", NULL, module_doc); 576 if (m == NULL) 577 return; 578 Py_INCREF(&Random_Type); 579 PyModule_AddObject(m, "Random", (PyObject *)&Random_Type); 580}