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         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 = PyInt_FromLong((long)(self->state[i]));
 323                 if (element == NULL)
 324                         goto Fail;
 325                 PyTuple_SET_ITEM(state, i, element);
 326         }
 327         element = PyInt_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         long element;
 343
 344         if (!PyTuple_Check(state)) {
 345                 PyErr_SetString(PyExc_TypeError,
 346                         "state vector must be a tuple");
 347                 return NULL;
 348         }
 349         if (PyTuple_Size(state) != N+1) {
 350                 PyErr_SetString(PyExc_ValueError,
 351                         "state vector is the wrong size");
 352                 return NULL;
 353         }
 354
 355         for (i=0; i<N ; i++) {
 356                 element = PyInt_AsLong(PyTuple_GET_ITEM(state, i));
 357                 if (element == -1 && PyErr_Occurred())
 358                         return NULL;
 359                 self->state[i] = (unsigned long)element;
 360         }
 361
 362         element = PyInt_AsLong(PyTuple_GET_ITEM(state, i));
 363         if (element == -1 && PyErr_Occurred())
 364                 return NULL;
 365         self->index = (int)element;
 366
 367         Py_INCREF(Py_None);
 368         return Py_None;
 369 }
 370
 371 /*
 372 Jumpahead should be a fast way advance the generator n-steps ahead, but
 373 lacking a formula for that, the next best is to use n and the existing
 374 state to create a new state far away from the original.
 375
 376 The generator uses constant spaced additive feedback, so shuffling the
 377 state elements ought to produce a state which would not be encountered
 378 (in the near term) by calls to random().  Shuffling is normally
 379 implemented by swapping the ith element with another element ranging
 380 from 0 to i inclusive.  That allows the element to have the possibility
 381 of not being moved.  Since the goal is to produce a new, different
 382 state, the swap element is ranged from 0 to i-1 inclusive.  This assures
 383 that each element gets moved at least once.
 384
 385 To make sure that consecutive calls to jumpahead(n) produce different
 386 states (even in the rare case of involutory shuffles), i+1 is added to
 387 each element at position i.  Successive calls are then guaranteed to
 388 have changing (growing) values as well as shuffled positions.
 389
 390 Finally, the self->index value is set to N so that the generator itself
 391 kicks in on the next call to random().  This assures that all results
 392 have been through the generator and do not just reflect alterations to
 393 the underlying state.
 394 */
 395
 396 static PyObject *
 397 random_jumpahead(RandomObject *self, PyObject *n)
 398 {
 399         long i, j;
 400         PyObject *iobj;
 401         PyObject *remobj;
 402         unsigned long *mt, tmp;
 403
 404         if (!PyInt_Check(n) && !PyLong_Check(n)) {
 405                 PyErr_Format(PyExc_TypeError, "jumpahead requires an "
 406                              "integer, not '%s'",
 407                              n->ob_type->tp_name);
 408                 return NULL;
 409         }
 410
 411         mt = self->state;
 412         for (i = N-1; i > 1; i--) {
 413                 iobj = PyInt_FromLong(i);
 414                 if (iobj == NULL)
 415                         return NULL;
 416                 remobj = PyNumber_Remainder(n, iobj);
 417                 Py_DECREF(iobj);
 418                 if (remobj == NULL)
 419                         return NULL;
 420                 j = PyInt_AsLong(remobj);
 421                 Py_DECREF(remobj);
 422                 if (j == -1L && PyErr_Occurred())
 423                         return NULL;
 424                 tmp = mt[i];
 425                 mt[i] = mt[j];
 426                 mt[j] = tmp;
 427         }
 428
 429         for (i = 0; i < N; i++)
 430                 mt[i] += i+1;
 431
 432         self->index = N;
 433         Py_INCREF(Py_None);
 434         return Py_None;
 435 }
 436
 437 static PyObject *
 438 random_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
 439 {
 440         RandomObject *self;
 441         PyObject *tmp;
 442
 443         self = (RandomObject *)type->tp_alloc(type, 0);
 444         if (self == NULL)
 445                 return NULL;
 446         tmp = random_seed(self, args);
 447         if (tmp == NULL) {
 448                 Py_DECREF(self);
 449                 return NULL;
 450         }
 451         Py_DECREF(tmp);
 452         return (PyObject *)self;
 453 }
 454
 455 static PyMethodDef random_methods[] = {
 456         {"random",      (PyCFunction)random_random,  METH_NOARGS,
 457                 PyDoc_STR("random() -> x in the interval [0, 1).")},
 458         {"seed",        (PyCFunction)random_seed,  METH_VARARGS,
 459                 PyDoc_STR("seed([n]) -> None.  Defaults to current time.")},
 460         {"getstate",    (PyCFunction)random_getstate,  METH_NOARGS,
 461                 PyDoc_STR("getstate() -> tuple containing the current state.")},
 462         {"setstate",      (PyCFunction)random_setstate,  METH_O,
 463                 PyDoc_STR("setstate(state) -> None.  Restores generator state.")},
 464         {"jumpahead",   (PyCFunction)random_jumpahead,  METH_O,
 465                 PyDoc_STR("jumpahead(int) -> None.  Create new state from "
 466                           "existing state and integer.")},
 467         {NULL,          NULL}           /* sentinel */
 468 };
 469
 470 PyDoc_STRVAR(random_doc,
 471 "Random() -> create a random number generator with its own internal state.");
 472
 473 static PyTypeObject Random_Type = {
 474         PyObject_HEAD_INIT(NULL)
 475         0,                              /*ob_size*/
 476         "_random.Random",               /*tp_name*/
 477         sizeof(RandomObject),           /*tp_basicsize*/
 478         0,                              /*tp_itemsize*/
 479         /* methods */
 480         0,                              /*tp_dealloc*/
 481         0,                              /*tp_print*/
 482         0,                              /*tp_getattr*/
 483         0,                              /*tp_setattr*/
 484         0,                              /*tp_compare*/
 485         0,                              /*tp_repr*/
 486         0,                              /*tp_as_number*/
 487         0,                              /*tp_as_sequence*/
 488         0,                              /*tp_as_mapping*/
 489         0,                              /*tp_hash*/
 490         0,                              /*tp_call*/
 491         0,                              /*tp_str*/
 492         PyObject_GenericGetAttr,        /*tp_getattro*/
 493         0,                              /*tp_setattro*/
 494         0,                              /*tp_as_buffer*/
 495         Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE,       /*tp_flags*/
 496         random_doc,                     /*tp_doc*/
 497         0,                              /*tp_traverse*/
 498         0,                              /*tp_clear*/
 499         0,                              /*tp_richcompare*/
 500         0,                              /*tp_weaklistoffset*/
 501         0,                              /*tp_iter*/
 502         0,                              /*tp_iternext*/
 503         random_methods,                 /*tp_methods*/
 504         0,                              /*tp_members*/
 505         0,                              /*tp_getset*/
 506         0,                              /*tp_base*/
 507         0,                              /*tp_dict*/
 508         0,                              /*tp_descr_get*/
 509         0,                              /*tp_descr_set*/
 510         0,                              /*tp_dictoffset*/
 511         0,                              /*tp_init*/
 512         0,                              /*tp_alloc*/
 513         random_new,                     /*tp_new*/
 514         _PyObject_Del,                  /*tp_free*/
 515         0,                              /*tp_is_gc*/
 516 };
 517
 518 PyDoc_STRVAR(module_doc,
 519 "Module implements the Mersenne Twister random number generator.");
 520
 521 PyMODINIT_FUNC
 522 init_random(void)
 523 {
 524         PyObject *m;
 525
 526         if (PyType_Ready(&Random_Type) < 0)
 527                 return;
 528         m = Py_InitModule3("_random", NULL, module_doc);
 529         Py_INCREF(&Random_Type);
 530         PyModule_AddObject(m, "Random", (PyObject *)&Random_Type);
 531 }