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)      ((v)->ob_type == &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
 217         PyObject *arg = NULL;
 218
 219         if (!PyArg_UnpackTuple(args, "seed", 0, 1, &arg))
 220                 return NULL;
 221
 222         if (arg == NULL || arg == Py_None) {
 223                 time_t now;
 224
 225                 time(&now);
 226                 init_genrand(self, (unsigned long)now);
 227                 Py_INCREF(Py_None);
 228                 return Py_None;
 229         }
 230         /* If the arg is an int or long, use its absolute value; else use
 231          * the absolute value of its hash code.
 232          */
 233         if (PyInt_Check(arg) || PyLong_Check(arg))
 234                 n = PyNumber_Absolute(arg);
 235         else {
 236                 long hash = PyObject_Hash(arg);
 237                 if (hash == -1)
 238                         goto Done;
 239                 n = PyLong_FromUnsignedLong((unsigned long)hash);
 240         }
 241         if (n == NULL)
 242                 goto Done;
 243
 244         /* Now split n into 32-bit chunks, from the right.  Each piece is
 245          * stored into key, which has a capacity of keymax chunks, of which
 246          * keyused are filled.  Alas, the repeated shifting makes this a
 247          * quadratic-time algorithm; we'd really like to use
 248          * _PyLong_AsByteArray here, but then we'd have to break into the
 249          * long representation to figure out how big an array was needed
 250          * in advance.
 251          */
 252         keymax = 8;     /* arbitrary; grows later if needed */
 253         keyused = 0;
 254         key = (unsigned long *)PyMem_Malloc(keymax * sizeof(*key));
 255         if (key == NULL)
 256                 goto Done;
 257
 258         masklower = PyLong_FromUnsignedLong(0xffffffffU);
 259         if (masklower == NULL)
 260                 goto Done;
 261         thirtytwo = PyInt_FromLong(32L);
 262         if (thirtytwo == NULL)
 263                 goto Done;
 264         while (PyObject_IsTrue(n)) {
 265                 PyObject *newn;
 266                 PyObject *pychunk;
 267                 unsigned long chunk;
 268
 269                 pychunk = PyNumber_And(n, masklower);
 270                 if (pychunk == NULL)
 271                         goto Done;
 272                 chunk = PyLong_AsUnsignedLong(pychunk);
 273                 Py_DECREF(pychunk);
 274                 if (chunk == (unsigned long)-1 && PyErr_Occurred())
 275                         goto Done;
 276                 newn = PyNumber_Rshift(n, thirtytwo);
 277                 if (newn == NULL)
 278                         goto Done;
 279                 Py_DECREF(n);
 280                 n = newn;
 281                 if (keyused >= keymax) {
 282                         unsigned long bigger = keymax << 1;
 283                         if ((bigger >> 1) != keymax) {
 284                                 PyErr_NoMemory();
 285                                 goto Done;
 286                         }
 287                         key = (unsigned long *)PyMem_Realloc(key,
 288                                                 bigger * sizeof(*key));
 289                         if (key == NULL)
 290                                 goto Done;
 291                         keymax = bigger;
 292                 }
 293                 assert(keyused < keymax);
 294                 key[keyused++] = chunk;
 295         }
 296
 297         if (keyused == 0)
 298                 key[keyused++] = 0UL;
 299         result = init_by_array(self, key, keyused);
 300 Done:
 301         Py_XDECREF(masklower);
 302         Py_XDECREF(thirtytwo);
 303         Py_XDECREF(n);
 304         PyMem_Free(key);
 305         return result;
 306 }
 307
 308 static PyObject *
 309 random_getstate(RandomObject *self)
 310 {
 311         PyObject *state;
 312         PyObject *element;
 313         int i;
 314
 315         state = PyTuple_New(N+1);
 316         if (state == NULL)
 317                 return NULL;
 318         for (i=0; i<N ; i++) {
 319                 element = PyInt_FromLong((long)(self->state[i]));
 320                 if (element == NULL)
 321                         goto Fail;
 322                 PyTuple_SET_ITEM(state, i, element);
 323         }
 324         element = PyInt_FromLong((long)(self->index));
 325         if (element == NULL)
 326                 goto Fail;
 327         PyTuple_SET_ITEM(state, i, element);
 328         return state;
 329
 330 Fail:
 331         Py_DECREF(state);
 332         return NULL;
 333 }
 334
 335 static PyObject *
 336 random_setstate(RandomObject *self, PyObject *state)
 337 {
 338         int i;
 339         long element;
 340
 341         if (!PyTuple_Check(state)) {
 342                 PyErr_SetString(PyExc_TypeError,
 343                         "state vector must be a tuple");
 344                 return NULL;
 345         }
 346         if (PyTuple_Size(state) != N+1) {
 347                 PyErr_SetString(PyExc_ValueError,
 348                         "state vector is the wrong size");
 349                 return NULL;
 350         }
 351
 352         for (i=0; i<N ; i++) {
 353                 element = PyInt_AsLong(PyTuple_GET_ITEM(state, i));
 354                 if (element == -1 && PyErr_Occurred())
 355                         return NULL;
 356                 self->state[i] = (unsigned long)element;
 357         }
 358
 359         element = PyInt_AsLong(PyTuple_GET_ITEM(state, i));
 360         if (element == -1 && PyErr_Occurred())
 361                 return NULL;
 362         self->index = (int)element;
 363
 364         Py_INCREF(Py_None);
 365         return Py_None;
 366 }
 367
 368 /*
 369 Jumpahead should be a fast way advance the generator n-steps ahead, but
 370 lacking a formula for that, the next best is to use n and the existing
 371 state to create a new state far away from the original.
 372
 373 The generator uses constant spaced additive feedback, so shuffling the
 374 state elements ought to produce a state which would not be encountered
 375 (in the near term) by calls to random().  Shuffling is normally
 376 implemented by swapping the ith element with another element ranging
 377 from 0 to i inclusive.  That allows the element to have the possibility
 378 of not being moved.  Since the goal is to produce a new, different
 379 state, the swap element is ranged from 0 to i-1 inclusive.  This assures
 380 that each element gets moved at least once.
 381
 382 To make sure that consecutive calls to jumpahead(n) produce different
 383 states (even in the rare case of involutory shuffles), i+1 is added to
 384 each element at position i.  Successive calls are then guaranteed to
 385 have changing (growing) values as well as shuffled positions.
 386
 387 Finally, the self->index value is set to N so that the generator itself
 388 kicks in on the next call to random().  This assures that all results
 389 have been through the generator and do not just reflect alterations to
 390 the underlying state.
 391 */
 392
 393 static PyObject *
 394 random_jumpahead(RandomObject *self, PyObject *n)
 395 {
 396         long i, j;
 397         PyObject *iobj;
 398         PyObject *remobj;
 399         unsigned long *mt, tmp;
 400
 401         if (!PyInt_Check(n) && !PyLong_Check(n)) {
 402                 PyErr_Format(PyExc_TypeError, "jumpahead requires an "
 403                              "integer, not '%s'",
 404                              n->ob_type->tp_name);
 405                 return NULL;
 406         }
 407
 408         mt = self->state;
 409         for (i = N-1; i > 1; i--) {
 410                 iobj = PyInt_FromLong(i);
 411                 if (iobj == NULL)
 412                         return NULL;
 413                 remobj = PyNumber_Remainder(n, iobj);
 414                 Py_DECREF(iobj);
 415                 if (remobj == NULL)
 416                         return NULL;
 417                 j = PyInt_AsLong(remobj);
 418                 Py_DECREF(remobj);
 419                 if (j == -1L && PyErr_Occurred())
 420                         return NULL;
 421                 tmp = mt[i];
 422                 mt[i] = mt[j];
 423                 mt[j] = tmp;
 424         }
 425
 426         for (i = 0; i < N; i++)
 427                 mt[i] += i+1;
 428
 429         self->index = N;
 430         Py_INCREF(Py_None);
 431         return Py_None;
 432 }
 433
 434 static PyObject *
 435 random_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
 436 {
 437         RandomObject *self;
 438         PyObject *tmp;
 439
 440         self = (RandomObject *)type->tp_alloc(type, 0);
 441         if (self == NULL)
 442                 return NULL;
 443         tmp = random_seed(self, args);
 444         if (tmp == NULL) {
 445                 Py_DECREF(self);
 446                 return NULL;
 447         }
 448         Py_DECREF(tmp);
 449         return (PyObject *)self;
 450 }
 451
 452 static PyMethodDef random_methods[] = {
 453         {"random",      (PyCFunction)random_random,  METH_NOARGS,
 454                 PyDoc_STR("random() -> x in the interval [0, 1).")},
 455         {"seed",        (PyCFunction)random_seed,  METH_VARARGS,
 456                 PyDoc_STR("seed([n]) -> None.  Defaults to current time.")},
 457         {"getstate",    (PyCFunction)random_getstate,  METH_NOARGS,
 458                 PyDoc_STR("getstate() -> tuple containing the current state.")},
 459         {"setstate",      (PyCFunction)random_setstate,  METH_O,
 460                 PyDoc_STR("setstate(state) -> None.  Restores generator state.")},
 461         {"jumpahead",   (PyCFunction)random_jumpahead,  METH_O,
 462                 PyDoc_STR("jumpahead(int) -> None.  Create new state from "
 463                           "existing state and integer.")},
 464         {NULL,          NULL}           /* sentinel */
 465 };
 466
 467 PyDoc_STRVAR(random_doc,
 468 "Random() -> create a random number generator with its own internal state.");
 469
 470 static PyTypeObject Random_Type = {
 471         PyObject_HEAD_INIT(NULL)
 472         0,                              /*ob_size*/
 473         "_random.Random",               /*tp_name*/
 474         sizeof(RandomObject),           /*tp_basicsize*/
 475         0,                              /*tp_itemsize*/
 476         /* methods */
 477         0,                              /*tp_dealloc*/
 478         0,                              /*tp_print*/
 479         0,                              /*tp_getattr*/
 480         0,                              /*tp_setattr*/
 481         0,                              /*tp_compare*/
 482         0,                              /*tp_repr*/
 483         0,                              /*tp_as_number*/
 484         0,                              /*tp_as_sequence*/
 485         0,                              /*tp_as_mapping*/
 486         0,                              /*tp_hash*/
 487         0,                              /*tp_call*/
 488         0,                              /*tp_str*/
 489         PyObject_GenericGetAttr,        /*tp_getattro*/
 490         0,                              /*tp_setattro*/
 491         0,                              /*tp_as_buffer*/
 492         Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE,       /*tp_flags*/
 493         random_doc,                     /*tp_doc*/
 494         0,                              /*tp_traverse*/
 495         0,                              /*tp_clear*/
 496         0,                              /*tp_richcompare*/
 497         0,                              /*tp_weaklistoffset*/
 498         0,                              /*tp_iter*/
 499         0,                              /*tp_iternext*/
 500         random_methods,                 /*tp_methods*/
 501         0,                              /*tp_members*/
 502         0,                              /*tp_getset*/
 503         0,                              /*tp_base*/
 504         0,                              /*tp_dict*/
 505         0,                              /*tp_descr_get*/
 506         0,                              /*tp_descr_set*/
 507         0,                              /*tp_dictoffset*/
 508         0,                              /*tp_init*/
 509         PyType_GenericAlloc,            /*tp_alloc*/
 510         random_new,                     /*tp_new*/
 511         _PyObject_Del,                  /*tp_free*/
 512         0,                              /*tp_is_gc*/
 513 };
 514
 515 PyDoc_STRVAR(module_doc,
 516 "Module implements the Mersenne Twister random number generator.");
 517
 518 PyMODINIT_FUNC
 519 init_random(void)
 520 {
 521         PyObject *m;
 522
 523         if (PyType_Ready(&Random_Type) < 0)
 524                 return;
 525         m = Py_InitModule3("_random", NULL, module_doc);
 526         Py_INCREF(&Random_Type);
 527         PyModule_AddObject(m, "Random", (PyObject *)&Random_Type);
 528 }