Random number generator in Python core

Modules/_randommodule.c

--- a/Modules/_randommodule.c	Fri Dec 30 21:26:08 2011 +0100
+++ b/Modules/_randommodule.c	Tue Jan 03 21:48:43 2012 +0100
@@ -1,85 +1,13 @@
 /* Random objects */
-
-/* ------------------------------------------------------------------
-   The code in this module was based on a download from:
-      http://www.math.keio.ac.jp/~matumoto/MT2002/emt19937ar.html
-
-   It was modified in 2002 by Raymond Hettinger as follows:
-
-    * the principal computational lines untouched.
-
-    * renamed genrand_res53() to random_random() and wrapped
-      in python calling/return code.
-
-    * genrand_int32() and the helper functions, init_genrand()
-      and init_by_array(), were declared static, wrapped in
-      Python calling/return code.  also, their global data
-      references were replaced with structure references.
-
-    * unused functions from the original were deleted.
-      new, original C python code was added to implement the
-      Random() interface.
-
-   The following are the verbatim comments from the original code:
-
-   A C-program for MT19937, with initialization improved 2002/1/26.
-   Coded by Takuji Nishimura and Makoto Matsumoto.
-
-   Before using, initialize the state by using init_genrand(seed)
-   or init_by_array(init_key, key_length).
-
-   Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
-   All rights reserved.
-
-   Redistribution and use in source and binary forms, with or without
-   modification, are permitted provided that the following conditions
-   are met:
-
-     1. Redistributions of source code must retain the above copyright
-    notice, this list of conditions and the following disclaimer.
-
-     2. Redistributions in binary form must reproduce the above copyright
-    notice, this list of conditions and the following disclaimer in the
-    documentation and/or other materials provided with the distribution.
-
-     3. The names of its contributors may not be used to endorse or promote
-    products derived from this software without specific prior written
-    permission.
-
-   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
-   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
-   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
-   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
-   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
-   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
-   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
-   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
-   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
-   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
-   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-
-
-   Any feedback is very welcome.
-   http://www.math.keio.ac.jp/matumoto/emt.html
-   email: matumoto@math.keio.ac.jp
-*/
-
-/* ---------------------------------------------------------------*/
-
 #include "Python.h"
+#include "pyrandom.h"
 #include <time.h>               /* for seeding to current time */
 
-/* Period parameters -- These are all magic.  Don't change. */
-#define N 624
-#define M 397
-#define MATRIX_A 0x9908b0dfUL   /* constant vector a */
-#define UPPER_MASK 0x80000000UL /* most significant w-r bits */
-#define LOWER_MASK 0x7fffffffUL /* least significant r bits */
+#define N _Py_MT_N
 
 typedef struct {
     PyObject_HEAD
-    unsigned long state[N];
-    int index;
+    _Py_MT_RandomState rs;
 } RandomObject;
 
 static PyTypeObject Random_Type;
@@ -87,123 +15,19 @@
 #define RandomObject_Check(v)      (Py_TYPE(v) == &Random_Type)
 
 
-/* Random methods */
-
-
-/* generates a random number on [0,0xffffffff]-interval */
-static unsigned long
-genrand_int32(RandomObject *self)
-{
-    unsigned long y;
-    static unsigned long mag01[2]={0x0UL, MATRIX_A};
-    /* mag01[x] = x * MATRIX_A  for x=0,1 */
-    unsigned long *mt;
-
-    mt = self->state;
-    if (self->index >= N) { /* generate N words at one time */
-        int kk;
-
-        for (kk=0;kk<N-M;kk++) {
-            y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
-            mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];
-        }
-        for (;kk<N-1;kk++) {
-            y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
-            mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
-        }
-        y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
-        mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
-
-        self->index = 0;
-    }
-
-    y = mt[self->index++];
-    y ^= (y >> 11);
-    y ^= (y << 7) & 0x9d2c5680UL;
-    y ^= (y << 15) & 0xefc60000UL;
-    y ^= (y >> 18);
-    return y;
-}
-
-/* random_random is the function named genrand_res53 in the original code;
- * generates a random number on [0,1) with 53-bit resolution; note that
- * 9007199254740992 == 2**53; I assume they're spelling "/2**53" as
- * multiply-by-reciprocal in the (likely vain) hope that the compiler will
- * optimize the division away at compile-time.  67108864 is 2**26.  In
- * effect, a contains 27 random bits shifted left 26, and b fills in the
- * lower 26 bits of the 53-bit numerator.
- * The orginal code credited Isaku Wada for this algorithm, 2002/01/09.
- */
-static PyObject *
-random_random(RandomObject *self)
-{
-    unsigned long a=genrand_int32(self)>>5, b=genrand_int32(self)>>6;
-    return PyFloat_FromDouble((a*67108864.0+b)*(1.0/9007199254740992.0));
-}
-
-/* initializes mt[N] with a seed */
-static void
-init_genrand(RandomObject *self, unsigned long s)
-{
-    int mti;
-    unsigned long *mt;
-
-    mt = self->state;
-    mt[0]= s & 0xffffffffUL;
-    for (mti=1; mti<N; mti++) {
-        mt[mti] =
-        (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
-        /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
-        /* In the previous versions, MSBs of the seed affect   */
-        /* only MSBs of the array mt[].                                */
-        /* 2002/01/09 modified by Makoto Matsumoto                     */
-        mt[mti] &= 0xffffffffUL;
-        /* for >32 bit machines */
-    }
-    self->index = mti;
-    return;
-}
-
-/* initialize by an array with array-length */
-/* init_key is the array for initializing keys */
-/* key_length is its length */
-static PyObject *
-init_by_array(RandomObject *self, unsigned long init_key[], unsigned long key_length)
-{
-    unsigned int i, j, k;       /* was signed in the original code. RDH 12/16/2002 */
-    unsigned long *mt;
-
-    mt = self->state;
-    init_genrand(self, 19650218UL);
-    i=1; j=0;
-    k = (N>key_length ? N : key_length);
-    for (; k; k--) {
-        mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL))
-                 + init_key[j] + j; /* non linear */
-        mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
-        i++; j++;
-        if (i>=N) { mt[0] = mt[N-1]; i=1; }
-        if (j>=key_length) j=0;
-    }
-    for (k=N-1; k; k--) {
-        mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL))
-                 - i; /* non linear */
-        mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
-        i++;
-        if (i>=N) { mt[0] = mt[N-1]; i=1; }
-    }
-
-    mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
-    Py_INCREF(Py_None);
-    return Py_None;
-}
-
 /*
  * The rest is Python-specific code, neither part of, nor derived from, the
  * Twister download.
  */
 
 static PyObject *
+random_random(RandomObject *self)
+{
+    return PyFloat_FromDouble(_Py_MT_GenRand_Res53(&self->rs));
+}
+
+
+static PyObject *
 random_seed(RandomObject *self, PyObject *args)
 {
     PyObject *result = NULL;            /* guilty until proved innocent */
@@ -224,7 +48,7 @@
         time_t now;
 
         time(&now);
-        init_genrand(self, (unsigned long)now);
+        _Py_MT_GenRand_Init(&self->rs, (unsigned long)now);
         Py_INCREF(Py_None);
         return Py_None;
     }
@@ -299,7 +123,9 @@
 
     if (keyused == 0)
         key[keyused++] = 0UL;
-    result = init_by_array(self, key, keyused);
+    _Py_MT_GenRand_InitArray(&self->rs, key, keyused);
+    Py_INCREF(Py_None);
+    result = Py_None;
 Done:
     Py_XDECREF(masklower);
     Py_XDECREF(thirtytwo);
@@ -319,12 +145,12 @@
     if (state == NULL)
         return NULL;
     for (i=0; i<N ; i++) {
-        element = PyLong_FromUnsignedLong(self->state[i]);
+        element = PyLong_FromUnsignedLong(self->rs.state[i]);
         if (element == NULL)
             goto Fail;
         PyTuple_SET_ITEM(state, i, element);
     }
-    element = PyLong_FromLong((long)(self->index));
+    element = PyLong_FromLong((long)(self->rs.index));
     if (element == NULL)
         goto Fail;
     PyTuple_SET_ITEM(state, i, element);
@@ -357,13 +183,13 @@
         element = PyLong_AsUnsignedLong(PyTuple_GET_ITEM(state, i));
         if (element == (unsigned long)-1 && PyErr_Occurred())
             return NULL;
-        self->state[i] = element & 0xffffffffUL; /* Make sure we get sane state */
+        self->rs.state[i] = element & 0xffffffffUL; /* Make sure we get sane state */
     }
 
     index = PyLong_AsLong(PyTuple_GET_ITEM(state, i));
     if (index == -1 && PyErr_Occurred())
         return NULL;
-    self->index = (int)index;
+    self->rs.index = (int)index;
 
     Py_INCREF(Py_None);
     return Py_None;
@@ -395,7 +221,7 @@
 
     /* Fill-out whole words, byte-by-byte to avoid endianness issues */
     for (i=0 ; i<bytes ; i+=4, k-=32) {
-        r = genrand_int32(self);
+        r = _Py_MT_GenRand_Int32(&self->rs);
         if (k < 32)
             r >>= (32 - k);
         bytearray[i+0] = (unsigned char)r;