diff options
Diffstat (limited to 'numpy')
| -rw-r--r-- | numpy/core/include/numpy/npy_common.h | 3 | ||||
| -rw-r--r-- | numpy/random/mtrand/mt_compat.h | 68 | ||||
| -rw-r--r-- | numpy/random/mtrand/mtrand.pyx | 390 | ||||
| -rw-r--r-- | numpy/random/mtrand/numpy.pxd | 16 | ||||
| -rw-r--r-- | numpy/random/mtrand/randomkit.c | 224 | ||||
| -rw-r--r-- | numpy/random/mtrand/randomkit.h | 41 |
6 files changed, 697 insertions, 45 deletions
diff --git a/numpy/core/include/numpy/npy_common.h b/numpy/core/include/numpy/npy_common.h index 47ef94c92..baf5549d9 100644 --- a/numpy/core/include/numpy/npy_common.h +++ b/numpy/core/include/numpy/npy_common.h @@ -7,6 +7,9 @@ #include <npy_config.h> #endif +/* need Python.h for npy_intp, npy_uintp */ +#include <Python.h> + /* * gcc does not unroll even with -O3 * use with care, unrolling on modern cpus rarely speeds things up 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); |