Modules/_randommodule.c

   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
  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 == -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;
 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
 438 static PyObject *
 439 random_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
 479 static PyObject *
 480 random_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
 500 static 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
 518 PyDoc_STRVAR(random_doc,
 519 "Random() -> create a random number generator with its own internal state.");
 520
 521 static 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
 565 PyDoc_STRVAR(module_doc,
 566 "Module implements the Mersenne Twister random number generator.");
 567
 568 PyMODINIT_FUNC
 569 init_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 }