ENH: Add dtype argument to random.randint.

Random ndarrays of the following types can now be generated:

* np.bool,
* np.int8, np.uint8,
* np.int16, np.uint16,
* np.int32, np.uint32,
* np.int64, np.uint64,
* np.int_ (long), np.intp

The specification is by precision rather than by C type. Hence, on some
platforms np.int64 may be a `long` instead of `long long` even if the
specified dtype is `long long` because the two may have the same
precision. The resulting type depends on which c type numpy uses for the
given precision. The byteorder specification is also ignored, the
generated arrays are always in native byte order.

The dtype of the result could be made more explicit if desired without
changing the user visible results.

5 files changed, 694 insertions, 45 deletions

diff --git a/numpy/random/mtrand/mt_compat.h b/numpy/random/mtrand/mt_compat.h
new file mode 100644
index 000000000..ab56a553c
--- /dev/null
+++ b/numpy/random/mtrand/mt_compat.h
@@ -0,0 +1,68 @@
+/*
+ * This is a convenience header file providing compatibility utilities
+ * for supporting Python 2 and Python 3 in the same code base.
+ *
+ * It can be removed when Python 2.6 is dropped as PyCapsule is available
+ * in both Python 3.1+ and Python 2.7.
+ */
+
+#ifndef _MT_COMPAT_H_
+#define _MT_COMPAT_H_
+
+#include <Python.h>
+#include <numpy/npy_common.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
+/*
+ * PyCObject functions adapted to PyCapsules.
+ *
+ * The main job here is to get rid of the improved error handling
+ * of PyCapsules. It's a shame...
+ */
+#if PY_VERSION_HEX >= 0x03000000
+
+static NPY_INLINE PyObject *
+NpyCapsule_FromVoidPtr(void *ptr, void (*dtor)(PyObject *))
+{
+    PyObject *ret = PyCapsule_New(ptr, NULL, dtor);
+    if (ret == NULL) {
+        PyErr_Clear();
+    }
+    return ret;
+}
+
+static NPY_INLINE void *
+NpyCapsule_AsVoidPtr(PyObject *obj)
+{
+    void *ret = PyCapsule_GetPointer(obj, NULL);
+    if (ret == NULL) {
+        PyErr_Clear();
+    }
+    return ret;
+}
+
+#else
+
+static NPY_INLINE PyObject *
+NpyCapsule_FromVoidPtr(void *ptr, void (*dtor)(void *))
+{
+    return PyCObject_FromVoidPtr(ptr, dtor);
+}
+
+static NPY_INLINE void *
+NpyCapsule_AsVoidPtr(PyObject *ptr)
+{
+    return PyCObject_AsVoidPtr(ptr);
+}
+
+#endif
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* _COMPAT_H_ */
diff --git a/numpy/random/mtrand/mtrand.pyx b/numpy/random/mtrand/mtrand.pyx
index 5120857d0..b83b2a588 100644
--- a/numpy/random/mtrand/mtrand.pyx
+++ b/numpy/random/mtrand/mtrand.pyx
@@ -67,6 +67,17 @@ cdef extern from "randomkit.h":
     rk_error rk_altfill(void *buffer, size_t size, int strong,
             rk_state *state) nogil
     double rk_gauss(rk_state *state) nogil
+    void rk_random_uint64(npy_uint64 off, npy_uint64 rng, npy_intp cnt,
+                          npy_uint64 *out, rk_state *state) nogil
+    void rk_random_uint32(npy_uint32 off, npy_uint32 rng, npy_intp cnt,
+                          npy_uint32 *out, rk_state *state) nogil
+    void rk_random_uint16(npy_uint16 off, npy_uint16 rng, npy_intp cnt,
+                          npy_uint16 *out, rk_state *state) nogil
+    void rk_random_uint8(npy_uint8 off, npy_uint8 rng, npy_intp cnt,
+                         npy_uint8 *out, rk_state *state) nogil
+    void rk_random_bool(npy_bool off, npy_bool rng, npy_intp cnt,
+                        npy_bool *out, rk_state *state) nogil
+
 
 cdef extern from "distributions.h":
     # do not need the GIL, but they do need a lock on the state !! */
@@ -131,6 +142,7 @@ cimport cython
 import numpy as np
 import operator
 import warnings
+
 try:
     from threading import Lock
 except ImportError:
@@ -574,6 +586,304 @@ def _shape_from_size(size, d):
            shape = tuple(size) + (d,)
     return shape
 
+
+# Set up dictionary of integer types and relevant functions.
+#
+# The dictionary is keyed by dtype(...).name and the values
+# are a tuple (low, high, function), where low and high are
+# the bounds of the largest half open interval `[low, high)`
+# and the function is the relevant function to call for
+# that precision.
+#
+# The functions are all the same except for changed types in
+# a few places. It would be easy to template them.
+
+def _rand_bool(low, high, size, rngstate):
+    """
+    _rand_bool(low, high, size, rngstate)
+
+    See `_rand_int32` for documentation, only the return type changes.
+
+    """
+    cdef npy_bool off, rng, buf
+    cdef npy_bool *out
+    cdef ndarray array "arrayObject"
+    cdef npy_intp cnt
+    cdef rk_state *state = <rk_state *>NpyCapsule_AsVoidPtr(rngstate)
+
+    rng = <npy_bool>(high - low)
+    off = <npy_bool>(low)
+    if size is None:
+        rk_random_bool(off, rng, 1, &buf, state)
+        return buf
+    else:
+        array = <ndarray>np.empty(size, np.bool_)
+        cnt = PyArray_SIZE(array)
+        out = <npy_bool *>PyArray_DATA(array)
+        with nogil:
+            rk_random_bool(off, rng, cnt, out, state)
+        return array
+
+
+def _rand_int8(low, high, size, rngstate):
+    """
+    _rand_int8(low, high, size, rngstate)
+
+    See `_rand_int32` for documentation, only the return type changes.
+
+    """
+    cdef npy_uint8 off, rng, buf
+    cdef npy_uint8 *out
+    cdef ndarray array "arrayObject"
+    cdef npy_intp cnt
+    cdef rk_state *state = <rk_state *>NpyCapsule_AsVoidPtr(rngstate)
+
+    rng = <npy_uint8>(high - low)
+    off = <npy_uint8>(<npy_int8>low)
+    if size is None:
+        rk_random_uint8(off, rng, 1, &buf, state)
+        return <npy_int8>buf
+    else:
+        array = <ndarray>np.empty(size, np.int8)
+        cnt = PyArray_SIZE(array)
+        out = <npy_uint8 *>PyArray_DATA(array)
+        with nogil:
+            rk_random_uint8(off, rng, cnt, out, state)
+        return array
+
+
+def _rand_int16(low, high, size, rngstate):
+    """
+    _rand_int16(low, high, size, rngstate)
+
+    See `_rand_int32` for documentation, only the return type changes.
+
+    """
+    cdef npy_uint16 off, rng, buf
+    cdef npy_uint16 *out
+    cdef ndarray array "arrayObject"
+    cdef npy_intp cnt
+    cdef rk_state *state = <rk_state *>NpyCapsule_AsVoidPtr(rngstate)
+
+    rng = <npy_uint16>(high - low)
+    off = <npy_uint16>(<npy_int16>low)
+    if size is None:
+        rk_random_uint16(off, rng, 1, &buf, state)
+        return <npy_int16>buf
+    else:
+        array = <ndarray>np.empty(size, np.int16)
+        cnt = PyArray_SIZE(array)
+        out = <npy_uint16 *>PyArray_DATA(array)
+        with nogil:
+            rk_random_uint16(off, rng, cnt, out, state)
+        return array
+
+
+def _rand_int32(low, high, size, rngstate):
+    """
+    _rand_int32(self, low, high, size, rngstate)
+
+    Return random np.int32 integers between `low` and `high`, inclusive.
+
+    Return random integers from the "discrete uniform" distribution in the
+    closed interval [`low`, `high`).  If `high` is None (the default),
+    then results are from [0, `low`). On entry the arguments are presumed
+    to have been validated for size and order for the np.int32 type.
+
+    Parameters
+    ----------
+    low : int
+        Lowest (signed) integer to be drawn from the distribution (unless
+        ``high=None``, in which case this parameter is the *highest* such
+        integer).
+    high : int
+        If provided, the largest (signed) integer to be drawn from the
+        distribution (see above for behavior if ``high=None``).
+    size : int or tuple of ints
+        Output shape.  If the given shape is, e.g., ``(m, n, k)``, then
+        ``m * n * k`` samples are drawn.  Default is None, in which case a
+        single value is returned.
+    rngstate : encapsulated pointer to rk_state
+        The specific type depends on the python version. In Python 2 it is
+        a PyCObject, in Python 3 a PyCapsule object.
+
+    Returns
+    -------
+    out : python scalar or ndarray of np.int32
+          `size`-shaped array of random integers from the appropriate
+          distribution, or a single such random int if `size` not provided.
+
+    """
+    cdef npy_uint32 off, rng, buf
+    cdef npy_uint32 *out
+    cdef ndarray array "arrayObject"
+    cdef npy_intp cnt
+    cdef rk_state *state = <rk_state *>NpyCapsule_AsVoidPtr(rngstate)
+
+    rng = <npy_uint32>(high - low)
+    off = <npy_uint32>(<npy_int32>low)
+    if size is None:
+        rk_random_uint32(off, rng, 1, &buf, state)
+        return <npy_int32>buf
+    else:
+        array = <ndarray>np.empty(size, np.int32)
+        cnt = PyArray_SIZE(array)
+        out = <npy_uint32 *>PyArray_DATA(array)
+        with nogil:
+            rk_random_uint32(off, rng, cnt, out, state)
+        return array
+
+
+def _rand_int64(low, high, size, rngstate):
+    """
+    _rand_int64(low, high, size, rngstate)
+
+    See `_rand_int32` for documentation, only the return type changes.
+
+    """
+    cdef npy_uint64 off, rng, buf
+    cdef npy_uint64 *out
+    cdef ndarray array "arrayObject"
+    cdef npy_intp cnt
+    cdef rk_state *state = <rk_state *>NpyCapsule_AsVoidPtr(rngstate)
+
+    rng = <npy_uint64>(high - low)
+    off = <npy_uint64>(<npy_int64>low)
+    if size is None:
+        rk_random_uint64(off, rng, 1, &buf, state)
+        return <npy_int64>buf
+    else:
+        array = <ndarray>np.empty(size, np.int64)
+        cnt = PyArray_SIZE(array)
+        out = <npy_uint64 *>PyArray_DATA(array)
+        with nogil:
+            rk_random_uint64(off, rng, cnt, out, state)
+        return array
+
+def _rand_uint8(low, high, size, rngstate):
+    """
+    _rand_uint8(low, high, size, rngstate)
+
+    See `_rand_int32` for documentation, only the return type changes.
+
+    """
+    cdef npy_uint8 off, rng, buf
+    cdef npy_uint8 *out
+    cdef ndarray array "arrayObject"
+    cdef npy_intp cnt
+    cdef rk_state *state = <rk_state *>NpyCapsule_AsVoidPtr(rngstate)
+
+    rng = <npy_uint8>(high - low)
+    off = <npy_uint8>(low)
+    if size is None:
+        rk_random_uint8(off, rng, 1, &buf, state)
+        return buf
+    else:
+        array = <ndarray>np.empty(size, np.uint8)
+        cnt = PyArray_SIZE(array)
+        out = <npy_uint8 *>PyArray_DATA(array)
+        with nogil:
+            rk_random_uint8(off, rng, cnt, out, state)
+        return array
+
+
+def _rand_uint16(low, high, size, rngstate):
+    """
+    _rand_uint16(low, high, size, rngstate)
+
+    See `_rand_int32` for documentation, only the return type changes.
+
+    """
+    cdef npy_uint16 off, rng, buf
+    cdef npy_uint16 *out
+    cdef ndarray array "arrayObject"
+    cdef npy_intp cnt
+    cdef rk_state *state = <rk_state *>NpyCapsule_AsVoidPtr(rngstate)
+
+    rng = <npy_uint16>(high - low)
+    off = <npy_uint16>(low)
+    if size is None:
+        rk_random_uint16(off, rng, 1, &buf, state)
+        return buf
+    else:
+        array = <ndarray>np.empty(size, np.uint16)
+        cnt = PyArray_SIZE(array)
+        out = <npy_uint16 *>PyArray_DATA(array)
+        with nogil:
+            rk_random_uint16(off, rng, cnt, out, state)
+        return array
+
+
+def _rand_uint32(low, high, size, rngstate):
+    """
+    _rand_uint32(self, low, high, size, rngstate)
+
+    See `_rand_int32` for documentation, only the return type changes.
+
+    """
+    cdef npy_uint32 off, rng, buf
+    cdef npy_uint32 *out
+    cdef ndarray array "arrayObject"
+    cdef npy_intp cnt
+    cdef rk_state *state = <rk_state *>NpyCapsule_AsVoidPtr(rngstate)
+
+    rng = <npy_uint32>(high - low)
+    off = <npy_uint32>(low)
+    if size is None:
+        rk_random_uint32(off, rng, 1, &buf, state)
+        return <npy_uint32>buf
+    else:
+        array = <ndarray>np.empty(size, np.uint32)
+        cnt = PyArray_SIZE(array)
+        out = <npy_uint32 *>PyArray_DATA(array)
+        with nogil:
+            rk_random_uint32(off, rng, cnt, out, state)
+        return array
+
+
+def _rand_uint64(low, high, size, rngstate):
+    """
+    _rand_uint64(low, high, size, rngstate)
+
+    See `_rand_int32` for documentation, only the return type changes.
+
+    """
+    cdef npy_uint64 off, rng, buf
+    cdef npy_uint64 *out
+    cdef ndarray array "arrayObject"
+    cdef npy_intp cnt
+    cdef rk_state *state = <rk_state *>NpyCapsule_AsVoidPtr(rngstate)
+
+    rng = <npy_uint64>(high - low)
+    off = <npy_uint64>(low)
+    if size is None:
+        rk_random_uint64(off, rng, 1, &buf, state)
+        return <npy_uint64>buf
+    else:
+        array = <ndarray>np.empty(size, np.uint64)
+        cnt = PyArray_SIZE(array)
+        out = <npy_uint64 *>PyArray_DATA(array)
+        with nogil:
+            rk_random_uint64(off, rng, cnt, out, state)
+        return array
+
+# Look up table for randint functions keyed by type name. The stored data
+# is a tuple (lbnd, ubnd, func), where lbnd is the smallest value for the
+# type, ubnd is one greater than the largest value, and func is the
+# function to call.
+_randint_type = {
+    'bool': (0, 2, _rand_bool),
+    'int8': (-2**7, 2**7, _rand_int8),
+    'int16': (-2**15, 2**15, _rand_int16),
+    'int32': (-2**31, 2**31, _rand_int32),
+    'int64': (-2**63, 2**63, _rand_int64),
+    'uint8': (0, 2**8, _rand_uint8),
+    'uint16': (0, 2**16, _rand_uint16),
+    'uint32': (0, 2**32, _rand_uint32),
+    'uint64': (0, 2**64, _rand_uint64)
+    }
+
+
 cdef class RandomState:
     """
     RandomState(seed=None)
@@ -618,11 +928,12 @@ cdef class RandomState:
     """
     cdef rk_state *internal_state
     cdef object lock
+    cdef object state_address
     poisson_lam_max = np.iinfo('l').max - np.sqrt(np.iinfo('l').max)*10
 
     def __init__(self, seed=None):
         self.internal_state = <rk_state*>PyMem_Malloc(sizeof(rk_state))
-
+        self.state_address = NpyCapsule_FromVoidPtr(self.internal_state, NULL)
         self.lock = Lock()
         self.seed(seed)
 
@@ -885,15 +1196,15 @@ cdef class RandomState:
         """
         return disc0_array(self.internal_state, rk_long, size, self.lock)
 
-    def randint(self, low, high=None, size=None):
+    def randint(self, low, high=None, size=None, dtype='l'):
         """
-        randint(low, high=None, size=None)
+        randint(low, high=None, size=None, dtype='l')
 
         Return random integers from `low` (inclusive) to `high` (exclusive).
 
-        Return random integers from the "discrete uniform" distribution in the
-        "half-open" interval [`low`, `high`). If `high` is None (the default),
-        then results are from [0, `low`).
+        Return random integers from the "discrete uniform" distribution of
+        the specified dtype in the "half-open" interval [`low`, `high`). If
+        `high` is None (the default), then results are from [0, `low`).
 
         Parameters
         ----------
@@ -908,6 +1219,13 @@ cdef class RandomState:
             Output shape.  If the given shape is, e.g., ``(m, n, k)``, then
             ``m * n * k`` samples are drawn.  Default is None, in which case a
             single value is returned.
+        dtype : dtype, optional
+            Desired dtype of the result. All dtypes are determined by their
+            name, i.e., 'int64', 'int`, etc, so byteorder is not available
+            and a specific precision may have different C types depending
+            on the platform. The default value is 'l' (C long).
+
+            .. versionadded:: 1.11.0
 
         Returns
         -------
@@ -936,14 +1254,24 @@ cdef class RandomState:
                [3, 2, 2, 0]])
 
         """
-        if high is not None and low >= high:
-            raise ValueError("low >= high")
-
         if high is None:
             high = low
             low = 0
 
-        return self.random_integers(low, high - 1, size)
+        key = np.dtype(dtype).name
+        if not key in _randint_type:
+            raise TypeError('Unsupported dtype "%s" for randint' % key)
+        lowbnd, highbnd, randfunc = _randint_type[key]
+
+        if low < lowbnd:
+            raise ValueError("low is out of bounds for %s" % (key,))
+        if high > highbnd:
+            raise ValueError("high is out of bounds for %s" % (key,))
+        if low >= high:
+            raise ValueError("low >= high")
+
+        with self.lock:
+            return randfunc(low, high - 1, size, self.state_address)
 
     def bytes(self, npy_intp length):
         """
@@ -1356,11 +1684,13 @@ cdef class RandomState:
         """
         random_integers(low, high=None, size=None)
 
-        Return random integers between `low` and `high`, inclusive.
+        Random integers of type np.int between `low` and `high`, inclusive.
 
-        Return random integers from the "discrete uniform" distribution in the
-        closed interval [`low`, `high`].  If `high` is None (the default),
-        then results are from [1, `low`].
+        Return random integers of type np.int from the "discrete uniform"
+        distribution in the closed interval [`low`, `high`].  If `high` is
+        None (the default), then results are from [1, `low`]. The np.int
+        type translates to the C long type used by Python 2 for "short"
+        integers and its precision is platform dependent.
 
         Parameters
         ----------
@@ -1426,37 +1756,13 @@ cdef class RandomState:
         >>> plt.show()
 
         """
-        if high is not None and low > high:
-            raise ValueError("low > high")
+        if high is None:
+            high = low
+            low = 1
 
-        cdef long lo, hi, rv
-        cdef unsigned long diff
-        cdef long *array_data
-        cdef ndarray array "arrayObject"
-        cdef npy_intp length
-        cdef npy_intp i
+        return self.randint(low, high + 1, size=size, dtype='l')
 
-        if high is None:
-            lo = 1
-            hi = low
-        else:
-            lo = low
-            hi = high
 
-        diff = <unsigned long>hi - <unsigned long>lo
-        if size is None:
-            with self.lock:
-                rv = lo + <long>rk_interval(diff, self. internal_state)
-            return rv
-        else:
-            array = <ndarray>np.empty(size, int)
-            length = PyArray_SIZE(array)
-            array_data = <long *>PyArray_DATA(array)
-            with self.lock, nogil:
-                for i from 0 <= i < length:
-                    rv = lo + <long>rk_interval(diff, self. internal_state)
-                    array_data[i] = rv
-            return array
 
     # Complicated, continuous distributions:
     def standard_normal(self, size=None):
diff --git a/numpy/random/mtrand/numpy.pxd b/numpy/random/mtrand/numpy.pxd
index c54f79c0a..488278d6c 100644
--- a/numpy/random/mtrand/numpy.pxd
+++ b/numpy/random/mtrand/numpy.pxd
@@ -2,6 +2,12 @@
 
 cdef extern from "numpy/npy_no_deprecated_api.h": pass
 
+cdef extern from "mt_compat.h":
+
+    object NpyCapsule_FromVoidPtr(void *ptr, void (*dtor)(object o))
+    void * NpyCapsule_AsVoidPtr(object o)
+    
+
 cdef extern from "numpy/arrayobject.h":
 
     cdef enum NPY_TYPES:
@@ -71,7 +77,17 @@ cdef extern from "numpy/arrayobject.h":
         double real
         double imag
 
+    ctypedef int npy_int
     ctypedef int npy_intp
+    ctypedef int npy_int64
+    ctypedef int npy_uint64
+    ctypedef int npy_int32
+    ctypedef int npy_uint32
+    ctypedef int npy_int16
+    ctypedef int npy_uint16
+    ctypedef int npy_int8
+    ctypedef int npy_uint8
+    ctypedef int npy_bool
 
     ctypedef extern class numpy.dtype [object PyArray_Descr]: pass
 
diff --git a/numpy/random/mtrand/randomkit.c b/numpy/random/mtrand/randomkit.c
index b18897e2c..3a95efeeb 100644
--- a/numpy/random/mtrand/randomkit.c
+++ b/numpy/random/mtrand/randomkit.c
@@ -70,6 +70,7 @@
 #include <errno.h>
 #include <limits.h>
 #include <math.h>
+#include <assert.h>
 
 #ifdef _WIN32
 /*
@@ -115,6 +116,10 @@
 #include <unistd.h>
 #endif
 
+/*
+ * Do not move this include. randomkit.h must be included
+ * after windows timeb.h is included.
+ */
 #include "randomkit.h"
 
 #ifndef RK_DEV_URANDOM
@@ -207,7 +212,11 @@ rk_randomseed(rk_state *state)
 #define UPPER_MASK 0x80000000UL
 #define LOWER_MASK 0x7fffffffUL
 
-/* Slightly optimised reference implementation of the Mersenne Twister */
+/*
+ * Slightly optimised reference implementation of the Mersenne Twister
+ * Note that regardless of the precision of long, only 32 bit random
+ * integers are produced
+ */
 unsigned long
 rk_random(rk_state *state)
 {
@@ -240,6 +249,219 @@ rk_random(rk_state *state)
     return y;
 }
 
+
+/*
+ * Returns an unsigned 64 bit random integer.
+ */
+NPY_INLINE static npy_uint64
+rk_uint64(rk_state *state)
+{
+    npy_uint64 upper = (npy_uint64)rk_random(state) << 32;
+    npy_uint64 lower = (npy_uint64)rk_random(state);
+    return upper | lower;
+}
+
+
+/*
+ * Returns an unsigned 32 bit random integer.
+ */
+NPY_INLINE static npy_uint32
+rk_uint32(rk_state *state)
+{
+    return (npy_uint32)rk_random(state);
+}
+
+
+/*
+ * Fills an array with cnt random npy_uint64 between off and off + rng
+ * inclusive. The numbers wrap if rng is sufficiently large.
+ */
+void
+rk_random_uint64(npy_uint64 off, npy_uint64 rng, npy_intp cnt,
+                 npy_uint64 *out, rk_state *state)
+{
+    npy_uint64 val, mask = rng;
+    npy_intp i;
+
+    if (rng == 0) {
+        for (i = 0; i < cnt; i++) {
+            out[i] = off;
+        }
+        return;
+    }
+
+    /* Smallest bit mask >= max */
+    mask |= mask >> 1;
+    mask |= mask >> 2;
+    mask |= mask >> 4;
+    mask |= mask >> 8;
+    mask |= mask >> 16;
+    mask |= mask >> 32;
+
+    for (i = 0; i < cnt; i++) {
+        if (rng <= 0xffffffffUL) {
+            while ((val = (rk_uint32(state) & mask)) > rng);
+        }
+        else {
+            while ((val = (rk_uint64(state) & mask)) > rng);
+        }
+        out[i] =  off + val;
+    }
+}
+
+
+/*
+ * Fills an array with cnt random npy_uint32 between off and off + rng
+ * inclusive. The numbers wrap if rng is sufficiently large.
+ */
+void
+rk_random_uint32(npy_uint32 off, npy_uint32 rng, npy_intp cnt,
+                 npy_uint32 *out, rk_state *state)
+{
+    npy_uint32 val, mask = rng;
+    npy_intp i;
+
+    if (rng == 0) {
+        for (i = 0; i < cnt; i++) {
+            out[i] = off;
+        }
+        return;
+    }
+
+    /* Smallest bit mask >= max */
+    mask |= mask >> 1;
+    mask |= mask >> 2;
+    mask |= mask >> 4;
+    mask |= mask >> 8;
+    mask |= mask >> 16;
+
+    for (i = 0; i < cnt; i++) {
+        while ((val = (rk_uint32(state) & mask)) > rng);
+        out[i] =  off + val;
+    }
+}
+
+
+/*
+ * Fills an array with cnt random npy_uint16 between off and off + rng
+ * inclusive. The numbers wrap if rng is sufficiently large.
+ */
+void
+rk_random_uint16(npy_uint16 off, npy_uint16 rng, npy_intp cnt,
+                 npy_uint16 *out, rk_state *state)
+{
+    npy_uint16 val, mask = rng;
+    npy_intp i;
+    npy_uint32 buf;
+    int bcnt = 0;
+
+    if (rng == 0) {
+        for (i = 0; i < cnt; i++) {
+            out[i] = off;
+        }
+        return;
+    }
+
+    /* Smallest bit mask >= max */
+    mask |= mask >> 1;
+    mask |= mask >> 2;
+    mask |= mask >> 4;
+    mask |= mask >> 8;
+
+    for (i = 0; i < cnt; i++) {
+        do {
+            if (!bcnt) {
+                buf = rk_uint32(state);
+                bcnt = 1;
+            }
+            else {
+                buf >>= 16;
+                bcnt--;
+            }
+            val = (npy_uint16)buf & mask;
+        } while (val > rng);
+        out[i] =  off + val;
+    }
+}
+
+
+/*
+ * Fills an array with cnt random npy_uint8 between off and off + rng
+ * inclusive. The numbers wrap if rng is sufficiently large.
+ */
+void
+rk_random_uint8(npy_uint8 off, npy_uint8 rng, npy_intp cnt,
+                npy_uint8 *out, rk_state *state)
+{
+    npy_uint8 val, mask = rng;
+    npy_intp i;
+    npy_uint32 buf;
+    int bcnt = 0;
+
+    if (rng == 0) {
+        for (i = 0; i < cnt; i++) {
+            out[i] = off;
+        }
+        return;
+    }
+
+    /* Smallest bit mask >= max */
+    mask |= mask >> 1;
+    mask |= mask >> 2;
+    mask |= mask >> 4;
+
+    for (i = 0; i < cnt; i++) {
+        do {
+            if (!bcnt) {
+                buf = rk_uint32(state);
+                bcnt = 3;
+            }
+            else {
+                buf >>= 8;
+                bcnt--;
+            }
+            val = (npy_uint8)buf & mask;
+        } while (val > rng);
+        out[i] =  off + val;
+    }
+}
+
+
+/*
+ * Fills an array with cnt random npy_bool between off and off + rng
+ * inclusive.
+ */
+void
+rk_random_bool(npy_bool off, npy_bool rng, npy_intp cnt,
+                npy_bool *out, rk_state *state)
+{
+    npy_intp i;
+    npy_uint32 buf;
+    int bcnt = 0;
+
+    if (rng == 0) {
+        for (i = 0; i < cnt; i++) {
+            out[i] = off;
+        }
+        return;
+    }
+
+    /* If we reach here rng and mask are one and off is zero */
+    assert(rng == 1 && off == 0);
+    for (i = 0; i < cnt; i++) {
+        if (!bcnt) {
+            buf = rk_uint32(state);
+            bcnt = 31;
+        }
+        else {
+            buf >>= 1;
+            bcnt--;
+        }
+        out[i] = (buf & 0x00000001) != 0;
+    }
+}
+
+
 long
 rk_long(rk_state *state)
 {
diff --git a/numpy/random/mtrand/randomkit.h b/numpy/random/mtrand/randomkit.h
index e049488ee..fcdd606a1 100644
--- a/numpy/random/mtrand/randomkit.h
+++ b/numpy/random/mtrand/randomkit.h
@@ -56,11 +56,13 @@
  *                 defaults to "/dev/urandom"
  */
 
-#include <stddef.h>
-
 #ifndef _RANDOMKIT_
 #define _RANDOMKIT_
 
+#include <stddef.h>
+#include <numpy/npy_common.h>
+
+
 #define RK_STATE_LEN 624
 
 typedef struct rk_state_
@@ -149,6 +151,41 @@ extern unsigned long rk_ulong(rk_state *state);
 extern unsigned long rk_interval(unsigned long max, rk_state *state);
 
 /*
+ * Fills an array with cnt random npy_uint64 between off and off + rng
+ * inclusive. The numbers wrap if rng is sufficiently large.
+ */
+extern void rk_random_uint64(npy_uint64 off, npy_uint64 rng, npy_intp cnt,
+                             npy_uint64 *out, rk_state *state);
+
+/*
+ * Fills an array with cnt random npy_uint32 between off and off + rng
+ * inclusive. The numbers wrap if rng is sufficiently large.
+ */
+extern void rk_random_uint32(npy_uint32 off, npy_uint32 rng, npy_intp cnt,
+                             npy_uint32 *out, rk_state *state);
+
+/*
+ * Fills an array with cnt random npy_uint16 between off and off + rng
+ * inclusive. The numbers wrap if rng is sufficiently large.
+ */
+extern void rk_random_uint16(npy_uint16 off, npy_uint16 rng, npy_intp cnt,
+                             npy_uint16 *out, rk_state *state);
+
+/*
+ * Fills an array with cnt random npy_uint8 between off and off + rng
+ * inclusive. The numbers wrap if rng is sufficiently large.
+ */
+extern void rk_random_uint8(npy_uint8 off, npy_uint8 rng, npy_intp cnt,
+                            npy_uint8 *out, rk_state *state);
+
+/*
+ * Fills an array with cnt random npy_bool between off and off + rng
+ * inclusive. It is assumed tha npy_bool as the same size as npy_uint8.
+ */
+extern void rk_random_bool(npy_bool off, npy_bool rng, npy_intp cnt,
+                           npy_bool *out, rk_state *state);
+
+/*
  * Returns a random double between 0.0 and 1.0, 1.0 excluded.
  */
 extern double rk_double(rk_state *state);