diff options
63 files changed, 976 insertions, 5079 deletions
diff --git a/LICENSES_bundled.txt b/LICENSES_bundled.txt index 03d22521a..ea349c7ee 100644 --- a/LICENSES_bundled.txt +++ b/LICENSES_bundled.txt @@ -23,5 +23,5 @@ License: BSD derived Name: dragon4 Files: numpy/core/src/multiarray/dragon4.c -License: One of a kind +License: MIT For license text, see numpy/core/src/multiarray/dragon4.c diff --git a/benchmarks/benchmarks/bench_random.py b/benchmarks/benchmarks/bench_random.py index 621b72d3a..3e47d35eb 100644 --- a/benchmarks/benchmarks/bench_random.py +++ b/benchmarks/benchmarks/bench_random.py @@ -92,7 +92,7 @@ nom_size = 100000 class RNG(Benchmark): param_names = ['rng'] - params = ['DSFMT', 'PCG64', 'PCG32', 'MT19937', 'Xoshiro256', + params = ['PCG64', 'PCG32', 'MT19937', 'Xoshiro256', 'Xoshiro512', 'Philox', 'ThreeFry', 'numpy'] def setup(self, bitgen): @@ -134,7 +134,7 @@ class Bounded(Benchmark): u32 = np.uint32 u64 = np.uint64 param_names = ['rng', 'dt_max'] - params = [['DSFMT', 'PCG64', 'PCG32', 'MT19937','Xoshiro256', + params = [['PCG64', 'PCG32', 'MT19937','Xoshiro256', 'Xoshiro512', 'Philox', 'ThreeFry', 'numpy'], [[u8, 95], [u8, 64], # Worst case for legacy diff --git a/doc/source/reference/arrays.dtypes.rst b/doc/source/reference/arrays.dtypes.rst index b55feb247..ab743a8ee 100644 --- a/doc/source/reference/arrays.dtypes.rst +++ b/doc/source/reference/arrays.dtypes.rst @@ -538,6 +538,7 @@ Attributes providing additional information: dtype.isnative dtype.descr dtype.alignment + dtype.base Methods diff --git a/doc/source/reference/c-api.types-and-structures.rst b/doc/source/reference/c-api.types-and-structures.rst index b72d9f902..a716b5a06 100644 --- a/doc/source/reference/c-api.types-and-structures.rst +++ b/doc/source/reference/c-api.types-and-structures.rst @@ -57,8 +57,8 @@ types are place holders that allow the array scalars to fit into a hierarchy of actual Python types. -PyArray_Type ------------- +PyArray_Type and PyArrayObject +------------------------------ .. c:var:: PyArray_Type @@ -74,7 +74,7 @@ PyArray_Type subclasses) will have this structure. For future compatibility, these structure members should normally be accessed using the provided macros. If you need a shorter name, then you can make use - of :c:type:`NPY_AO` which is defined to be equivalent to + of :c:type:`NPY_AO` (deprecated) which is defined to be equivalent to :c:type:`PyArrayObject`. .. code-block:: c @@ -91,7 +91,7 @@ PyArray_Type PyObject *weakreflist; } PyArrayObject; -.. c:macro: PyArrayObject.PyObject_HEAD +.. c:macro:: PyArrayObject.PyObject_HEAD This is needed by all Python objects. It consists of (at least) a reference count member ( ``ob_refcnt`` ) and a pointer to the @@ -130,14 +130,16 @@ PyArray_Type .. c:member:: PyObject *PyArrayObject.base This member is used to hold a pointer to another Python object that - is related to this array. There are two use cases: 1) If this array - does not own its own memory, then base points to the Python object - that owns it (perhaps another array object), 2) If this array has - the (deprecated) :c:data:`NPY_ARRAY_UPDATEIFCOPY` or - :c:data:NPY_ARRAY_WRITEBACKIFCOPY`: flag set, then this array is - a working copy of a "misbehaved" array. When - ``PyArray_ResolveWritebackIfCopy`` is called, the array pointed to by base - will be updated with the contents of this array. + is related to this array. There are two use cases: + + - If this array does not own its own memory, then base points to the + Python object that owns it (perhaps another array object) + - If this array has the (deprecated) :c:data:`NPY_ARRAY_UPDATEIFCOPY` or + :c:data:`NPY_ARRAY_WRITEBACKIFCOPY` flag set, then this array is a working + copy of a "misbehaved" array. + + When ``PyArray_ResolveWritebackIfCopy`` is called, the array pointed to + by base will be updated with the contents of this array. .. c:member:: PyArray_Descr *PyArrayObject.descr @@ -163,8 +165,8 @@ PyArray_Type weakref module). -PyArrayDescr_Type ------------------ +PyArrayDescr_Type and PyArray_Descr +----------------------------------- .. c:var:: PyArrayDescr_Type @@ -253,11 +255,13 @@ PyArrayDescr_Type .. c:var:: NPY_ITEM_REFCOUNT - .. c:var:: NPY_ITEM_HASOBJECT - Indicates that items of this data-type must be reference counted (using :c:func:`Py_INCREF` and :c:func:`Py_DECREF` ). + .. c:var:: NPY_ITEM_HASOBJECT + + Same as :c:data:`NPY_ITEM_REFCOUNT`. + .. c:var:: NPY_LIST_PICKLE Indicates arrays of this data-type must be converted to a list @@ -676,25 +680,28 @@ PyArrayDescr_Type The :c:data:`PyArray_Type` typeobject implements many of the features of -Python objects including the tp_as_number, tp_as_sequence, -tp_as_mapping, and tp_as_buffer interfaces. The rich comparison -(tp_richcompare) is also used along with new-style attribute lookup -for methods (tp_methods) and properties (tp_getset). The -:c:data:`PyArray_Type` can also be sub-typed. +:c:type:`Python objects <PyTypeObject>` including the :c:member:`tp_as_number +<PyTypeObject.tp_as_number>`, :c:member:`tp_as_sequence +<PyTypeObject.tp_as_sequence>`, :c:member:`tp_as_mapping +<PyTypeObject.tp_as_mapping>`, and :c:member:`tp_as_buffer +<PyTypeObject.tp_as_buffer>` interfaces. The :c:type:`rich comparison +<richcmpfunc>`) is also used along with new-style attribute lookup for +member (:c:member:`tp_members <PyTypeObject.tp_members>`) and properties +(:c:member:`tp_getset <PyTypeObject.tp_getset>`). +The :c:data:`PyArray_Type` can also be sub-typed. .. tip:: - The tp_as_number methods use a generic approach to call whatever - function has been registered for handling the operation. The - function PyNumeric_SetOps(..) can be used to register functions to - handle particular mathematical operations (for all arrays). When - the umath module is imported, it sets the numeric operations for - all arrays to the corresponding ufuncs. The tp_str and tp_repr - methods can also be altered using PyString_SetStringFunction(...). + The ``tp_as_number`` methods use a generic approach to call whatever + function has been registered for handling the operation. When the + ``_multiarray_umath module`` is imported, it sets the numeric operations + for all arrays to the corresponding ufuncs. This choice can be changed with + :c:func:`PyUFunc_ReplaceLoopBySignature` The ``tp_str`` and ``tp_repr`` + methods can also be altered using :c:func:`PyArray_SetStringFunction`. -PyUFunc_Type ------------- +PyUFunc_Type and PyUFuncObject +------------------------------ .. c:var:: PyUFunc_Type @@ -786,8 +793,8 @@ PyUFunc_Type the identity for this operation. It is only used for a reduce-like call on an empty array. - .. c:member:: void PyUFuncObject.functions(char** args, npy_intp* dims, - npy_intp* steps, void* extradata) + .. c:member:: void PyUFuncObject.functions( \ + char** args, npy_intp* dims, npy_intp* steps, void* extradata) An array of function pointers --- one for each data type supported by the ufunc. This is the vector loop that is called @@ -932,8 +939,8 @@ PyUFunc_Type - :c:data:`UFUNC_CORE_DIM_SIZE_INFERRED` if the dim size will be determined from the operands and not from a :ref:`frozen <frozen>` signature -PyArrayIter_Type ----------------- +PyArrayIter_Type and PyArrayIterObject +-------------------------------------- .. c:var:: PyArrayIter_Type @@ -1042,8 +1049,8 @@ with it through the use of the macros :c:func:`PyArray_ITER_NEXT` (it), :c:type:`PyArrayIterObject *`. -PyArrayMultiIter_Type ---------------------- +PyArrayMultiIter_Type and PyArrayMultiIterObject +------------------------------------------------ .. c:var:: PyArrayMultiIter_Type @@ -1104,8 +1111,8 @@ PyArrayMultiIter_Type arrays to be broadcast together. On return, the iterators are adjusted for broadcasting. -PyArrayNeighborhoodIter_Type ----------------------------- +PyArrayNeighborhoodIter_Type and PyArrayNeighborhoodIterObject +-------------------------------------------------------------- .. c:var:: PyArrayNeighborhoodIter_Type @@ -1118,8 +1125,33 @@ PyArrayNeighborhoodIter_Type :c:data:`PyArrayNeighborhoodIter_Type` is the :c:type:`PyArrayNeighborhoodIterObject`. -PyArrayFlags_Type ------------------ + .. code-block:: c + + typedef struct { + PyObject_HEAD + int nd_m1; + npy_intp index, size; + npy_intp coordinates[NPY_MAXDIMS] + npy_intp dims_m1[NPY_MAXDIMS]; + npy_intp strides[NPY_MAXDIMS]; + npy_intp backstrides[NPY_MAXDIMS]; + npy_intp factors[NPY_MAXDIMS]; + PyArrayObject *ao; + char *dataptr; + npy_bool contiguous; + npy_intp bounds[NPY_MAXDIMS][2]; + npy_intp limits[NPY_MAXDIMS][2]; + npy_intp limits_sizes[NPY_MAXDIMS]; + npy_iter_get_dataptr_t translate; + npy_intp nd; + npy_intp dimensions[NPY_MAXDIMS]; + PyArrayIterObject* _internal_iter; + char* constant; + int mode; + } PyArrayNeighborhoodIterObject; + +PyArrayFlags_Type and PyArrayFlagsObject +---------------------------------------- .. c:var:: PyArrayFlags_Type @@ -1129,6 +1161,16 @@ PyArrayFlags_Type attributes or by accessing them as if the object were a dictionary with the flag names as entries. +.. c:type:: PyArrayFlagsObject + + .. code-block:: c + + typedef struct PyArrayFlagsObject { + PyObject_HEAD + PyObject *arr; + int flags; + } PyArrayFlagsObject; + ScalarArrayTypes ---------------- diff --git a/doc/source/reference/random/bit_generators/dsfmt.rst b/doc/source/reference/random/bit_generators/dsfmt.rst deleted file mode 100644 index e7c6dbb31..000000000 --- a/doc/source/reference/random/bit_generators/dsfmt.rst +++ /dev/null @@ -1,36 +0,0 @@ -Double SIMD Mersenne Twister (dSFMT) ------------------------------------- - -.. module:: numpy.random.dsfmt - -.. currentmodule:: numpy.random.dsfmt - - -.. autoclass:: DSFMT - :exclude-members: - -Seeding and State -================= - -.. autosummary:: - :toctree: generated/ - - ~DSFMT.seed - ~DSFMT.state - -Parallel generation -=================== -.. autosummary:: - :toctree: generated/ - - ~DSFMT.jumped - -Extending -========= -.. autosummary:: - :toctree: generated/ - - ~DSFMT.cffi - ~DSFMT.ctypes - - diff --git a/doc/source/reference/random/bit_generators/index.rst b/doc/source/reference/random/bit_generators/index.rst index 3a9294bfb..4d3d39ae2 100644 --- a/doc/source/reference/random/bit_generators/index.rst +++ b/doc/source/reference/random/bit_generators/index.rst @@ -18,7 +18,6 @@ Stable RNGs .. toctree:: :maxdepth: 1 - DSFMT <dsfmt> MT19937 <mt19937> PCG32 <pcg32> PCG64 <pcg64> diff --git a/doc/source/reference/random/index.rst b/doc/source/reference/random/index.rst index 3159f0e1c..42956590a 100644 --- a/doc/source/reference/random/index.rst +++ b/doc/source/reference/random/index.rst @@ -153,10 +153,6 @@ The included BitGenerators are: * MT19937 - The standard Python BitGenerator. Produces identical results to Python using the same seed/state. Adds a `~mt19937.MT19937.jumped` function that returns a new generator with state as-if ``2**128`` draws have been made. -* dSFMT - SSE2 enabled versions of the MT19937 generator. Theoretically - the same, but with a different state and so it is not possible to produce a - sequence identical to MT19937. Supports ``jumped`` and so can - be used in parallel applications. See the `dSFMT authors' page`_. * Xorshiro256** and Xorshiro512** - The most recently introduced XOR, shift, and rotate generator. Supports ``jumped`` and so can be used in parallel applications. See the documentation for @@ -167,7 +163,6 @@ The included BitGenerators are: arbitrary number of steps or generating independent streams. See the `Random123`_ page for more details about this class of bit generators. -.. _`dSFMT authors' page`: http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/ .. _`PCG author's page`: http://www.pcg-random.org/ .. _`xorshift, xoroshiro and xoshiro authors' page`: http://xoroshiro.di.unimi.it/ .. _`Random123`: https://www.deshawresearch.com/resources_random123.html diff --git a/doc/source/reference/random/parallel.rst b/doc/source/reference/random/parallel.rst index ffbaea62b..6c495cc29 100644 --- a/doc/source/reference/random/parallel.rst +++ b/doc/source/reference/random/parallel.rst @@ -64,8 +64,6 @@ are listed below. +-----------------+-------------------------+-------------------------+-------------------------+ | BitGenerator | Period | Jump Size | Bits | +=================+=========================+=========================+=========================+ -| DSFMT | :math:`2^{19937}` | :math:`2^{128}` | 53 | -+-----------------+-------------------------+-------------------------+-------------------------+ | MT19937 | :math:`2^{19937}` | :math:`2^{128}` | 32 | +-----------------+-------------------------+-------------------------+-------------------------+ | PCG64 | :math:`2^{128}` | :math:`2^{64}` | 64 | diff --git a/doc/source/reference/random/performance.py b/doc/source/reference/random/performance.py index 54165226e..f1dc50c9d 100644 --- a/doc/source/reference/random/performance.py +++ b/doc/source/reference/random/performance.py @@ -4,10 +4,10 @@ from timeit import repeat import pandas as pd import numpy as np -from numpy.random import MT19937, DSFMT, ThreeFry, PCG64, Philox, \ +from numpy.random import MT19937, ThreeFry, PCG64, Philox, \ Xoshiro256, Xoshiro512 -PRNGS = [DSFMT, MT19937, PCG64, Philox, ThreeFry, Xoshiro256, Xoshiro512] +PRNGS = [MT19937, PCG64, Philox, ThreeFry, Xoshiro256, Xoshiro512] funcs = OrderedDict() integers = 'integers(0, 2**{bits},size=1000000, dtype="uint{bits}")' diff --git a/doc/source/reference/random/performance.rst b/doc/source/reference/random/performance.rst index 07867ee07..014db19a3 100644 --- a/doc/source/reference/random/performance.rst +++ b/doc/source/reference/random/performance.rst @@ -23,8 +23,7 @@ specific distribution. The original :class:`~mt19937.MT19937` generator is much slower since it requires 2 32-bit values to equal the output of the faster generators. -Integer performance has a similar ordering although `dSFMT` is slower since -it generates 53-bit floating point values rather than integer values. +Integer performance has a similar ordering. The pattern is similar for other, more complex generators. The normal performance of the legacy :class:`~mtrand.RandomState` generator is much @@ -120,9 +119,8 @@ Normal 115.3 100 135.6 60.3 93.6 1 ~~~~~~~~~~~~~~ The performance of 64-bit generators on 32-bit Windows is much lower than on 64-bit -operating systems due to register width. DSFMT uses SSE2 when available, and so is less -affected by the size of the operating system's register. MT19937, the generator that has been -in NumPy since 2005, operates on 32-bit integers and so is close to DSFMT. +operating systems due to register width. MT19937, the generator that has been +in NumPy since 2005, operates on 32-bit integers. =================== ======= ========= ======= ======== ========== ============ Distribution DSFMT MT19937 PCG64 Philox ThreeFry Xoshiro256 diff --git a/numpy/core/_add_newdocs.py b/numpy/core/_add_newdocs.py index 0b4bf3342..fafe29482 100644 --- a/numpy/core/_add_newdocs.py +++ b/numpy/core/_add_newdocs.py @@ -936,11 +936,11 @@ add_newdoc('numpy.core.multiarray', 'empty', -------- >>> np.empty([2, 2]) array([[ -9.74499359e+001, 6.69583040e-309], - [ 2.13182611e-314, 3.06959433e-309]]) #random + [ 2.13182611e-314, 3.06959433e-309]]) #uninitialized >>> np.empty([2, 2], dtype=int) array([[-1073741821, -1067949133], - [ 496041986, 19249760]]) #random + [ 496041986, 19249760]]) #uninitialized """) @@ -5390,6 +5390,17 @@ add_newdoc('numpy.core.multiarray', 'dtype', ('alignment', More information is available in the C-API section of the manual. + Examples + -------- + + >>> x = np.dtype('i4') + >>> x.alignment + 4 + + >>> x = np.dtype(float) + >>> x.alignment + 8 + """)) add_newdoc('numpy.core.multiarray', 'dtype', ('byteorder', @@ -5436,7 +5447,16 @@ add_newdoc('numpy.core.multiarray', 'dtype', ('byteorder', """)) add_newdoc('numpy.core.multiarray', 'dtype', ('char', - """A unique character code for each of the 21 different built-in types.""")) + """A unique character code for each of the 21 different built-in types. + + Examples + -------- + + >>> x = np.dtype(float) + >>> x.char + 'd' + + """)) add_newdoc('numpy.core.multiarray', 'dtype', ('descr', """ @@ -5447,6 +5467,18 @@ add_newdoc('numpy.core.multiarray', 'dtype', ('descr', Warning: This attribute exists specifically for `__array_interface__`, and is not a datatype description compatible with `np.dtype`. + + Examples + -------- + + >>> x = np.dtype(float) + >>> x.descr + [('', '<f8')] + + >>> dt = np.dtype([('name', np.str_, 16), ('grades', np.float64, (2,))]) + >>> dt.descr + [('name', '<U16'), ('grades', '<f8', (2,))] + """)) add_newdoc('numpy.core.multiarray', 'dtype', ('fields', @@ -5487,6 +5519,18 @@ add_newdoc('numpy.core.multiarray', 'dtype', ('flags', of these flags is in C-API documentation; they are largely useful for user-defined data-types. + The following example demonstrates that operations on this particular + dtype requires Python C-API. + + Examples + -------- + + >>> x = np.dtype([('a', np.int32, 8), ('b', np.float64, 6)]) + >>> x.flags + 16 + >>> np.core.multiarray.NEEDS_PYAPI + 16 + """)) add_newdoc('numpy.core.multiarray', 'dtype', ('hasobject', @@ -5544,6 +5588,7 @@ add_newdoc('numpy.core.multiarray', 'dtype', ('isalignedstruct', field alignment. This flag is sticky, so when combining multiple structs together, it is preserved and produces new dtypes which are also aligned. + """)) add_newdoc('numpy.core.multiarray', 'dtype', ('itemsize', @@ -5553,6 +5598,19 @@ add_newdoc('numpy.core.multiarray', 'dtype', ('itemsize', For 18 of the 21 types this number is fixed by the data-type. For the flexible data-types, this number can be anything. + Examples + -------- + + >>> arr = np.array([[1, 2], [3, 4]]) + >>> arr.dtype + dtype('int64') + >>> arr.itemsize + 8 + + >>> dt = np.dtype([('name', np.str_, 16), ('grades', np.float64, (2,))]) + >>> dt.itemsize + 80 + """)) add_newdoc('numpy.core.multiarray', 'dtype', ('kind', @@ -5573,6 +5631,19 @@ add_newdoc('numpy.core.multiarray', 'dtype', ('kind', V void = ====================== + Examples + -------- + + >>> dt = np.dtype('i4') + >>> dt.kind + 'i' + >>> dt = np.dtype('f8') + >>> dt.kind + 'f' + >>> dt = np.dtype([('field1', 'f8')]) + >>> dt.kind + 'V' + """)) add_newdoc('numpy.core.multiarray', 'dtype', ('name', @@ -5581,6 +5652,16 @@ add_newdoc('numpy.core.multiarray', 'dtype', ('name', Un-sized flexible data-type objects do not have this attribute. + Examples + -------- + + >>> x = np.dtype(float) + >>> x.name + 'float64' + >>> x = np.dtype([('a', np.int32, 8), ('b', np.float64, 6)]) + >>> x.name + 'void640' + """)) add_newdoc('numpy.core.multiarray', 'dtype', ('names', @@ -5604,6 +5685,17 @@ add_newdoc('numpy.core.multiarray', 'dtype', ('num', These are roughly ordered from least-to-most precision. + Examples + -------- + + >>> dt = np.dtype(str) + >>> dt.num + 19 + + >>> dt = np.dtype(float) + >>> dt.num + 12 + """)) add_newdoc('numpy.core.multiarray', 'dtype', ('shape', @@ -5611,6 +5703,17 @@ add_newdoc('numpy.core.multiarray', 'dtype', ('shape', Shape tuple of the sub-array if this data type describes a sub-array, and ``()`` otherwise. + Examples + -------- + + >>> dt = np.dtype(('i4', 4)) + >>> dt.shape + (4,) + + >>> dt = np.dtype(('i4', (2, 3))) + >>> dt.shape + (2, 3) + """)) add_newdoc('numpy.core.multiarray', 'dtype', ('ndim', @@ -5620,6 +5723,20 @@ add_newdoc('numpy.core.multiarray', 'dtype', ('ndim', .. versionadded:: 1.13.0 + Examples + -------- + >>> x = np.dtype(float) + >>> x.ndim + 0 + + >>> x = np.dtype((float, 8)) + >>> x.ndim + 1 + + >>> x = np.dtype(('i4', (3, 4))) + >>> x.ndim + 2 + """)) add_newdoc('numpy.core.multiarray', 'dtype', ('str', diff --git a/numpy/core/_internal.py b/numpy/core/_internal.py index ab5a64a1a..c70718cb6 100644 --- a/numpy/core/_internal.py +++ b/numpy/core/_internal.py @@ -146,7 +146,7 @@ def _reconstruct(subtype, shape, dtype): # format_re was originally from numarray by J. Todd Miller format_re = re.compile(br'(?P<order1>[<>|=]?)' - br'(?P<repeats> *[(]?[ ,0-9L]*[)]? *)' + br'(?P<repeats> *[(]?[ ,0-9]*[)]? *)' br'(?P<order2>[<>|=]?)' br'(?P<dtype>[A-Za-z0-9.?]*(?:\[[a-zA-Z0-9,.]+\])?)') sep_re = re.compile(br'\s*,\s*') diff --git a/numpy/core/code_generators/ufunc_docstrings.py b/numpy/core/code_generators/ufunc_docstrings.py index 2a5646bd7..fb418aadc 100644 --- a/numpy/core/code_generators/ufunc_docstrings.py +++ b/numpy/core/code_generators/ufunc_docstrings.py @@ -2601,8 +2601,8 @@ add_newdoc('numpy.core.umath', 'matmul', For other keyword-only arguments, see the :ref:`ufunc docs <ufuncs.kwargs>`. - ..versionadded:: 1.16 - Now handles ufunc kwargs + .. versionadded:: 1.16 + Now handles ufunc kwargs Returns ------- diff --git a/numpy/core/multiarray.py b/numpy/core/multiarray.py index 118562c89..4f2c5b78e 100644 --- a/numpy/core/multiarray.py +++ b/numpy/core/multiarray.py @@ -128,11 +128,11 @@ def empty_like(prototype, dtype=None, order=None, subok=None, shape=None): -------- >>> a = ([1,2,3], [4,5,6]) # a is array-like >>> np.empty_like(a) - array([[-1073741821, -1073741821, 3], # random + array([[-1073741821, -1073741821, 3], # uninitialized [ 0, 0, -1073741821]]) >>> a = np.array([[1., 2., 3.],[4.,5.,6.]]) >>> np.empty_like(a) - array([[ -2.00000715e+000, 1.48219694e-323, -2.00000572e+000], # random + array([[ -2.00000715e+000, 1.48219694e-323, -2.00000572e+000], # uninitialized [ 4.38791518e-305, -2.00000715e+000, 4.17269252e-309]]) """ diff --git a/numpy/core/src/multiarray/convert_datatype.c b/numpy/core/src/multiarray/convert_datatype.c index 8e24d6361..025c66013 100644 --- a/numpy/core/src/multiarray/convert_datatype.c +++ b/numpy/core/src/multiarray/convert_datatype.c @@ -1767,46 +1767,22 @@ dtype_kind_to_simplified_ordering(char kind) } } -/*NUMPY_API - * Produces the result type of a bunch of inputs, using the UFunc - * type promotion rules. Use this function when you have a set of - * input arrays, and need to determine an output array dtype. - * - * If all the inputs are scalars (have 0 dimensions) or the maximum "kind" - * of the scalars is greater than the maximum "kind" of the arrays, does - * a regular type promotion. - * - * Otherwise, does a type promotion on the MinScalarType - * of all the inputs. Data types passed directly are treated as array - * types. - * + +/* + * Determine if there is a mix of scalars and arrays/dtypes. + * If this is the case, the scalars should be handled as the minimum type + * capable of holding the value when the maximum "category" of the scalars + * surpasses the maximum "category" of the arrays/dtypes. + * If the scalars are of a lower or same category as the arrays, they may be + * demoted to a lower type within their category (the lowest type they can + * be cast to safely according to scalar casting rules). */ -NPY_NO_EXPORT PyArray_Descr * -PyArray_ResultType(npy_intp narrs, PyArrayObject **arr, - npy_intp ndtypes, PyArray_Descr **dtypes) +NPY_NO_EXPORT int +should_use_min_scalar(npy_intp narrs, PyArrayObject **arr, + npy_intp ndtypes, PyArray_Descr **dtypes) { - npy_intp i; - int use_min_scalar; - - /* If there's just one type, pass it through */ - if (narrs + ndtypes == 1) { - PyArray_Descr *ret = NULL; - if (narrs == 1) { - ret = PyArray_DESCR(arr[0]); - } - else { - ret = dtypes[0]; - } - Py_INCREF(ret); - return ret; - } + int use_min_scalar = 0; - /* - * Determine if there are any scalars, and if so, whether - * the maximum "kind" of the scalars surpasses the maximum - * "kind" of the arrays - */ - use_min_scalar = 0; if (narrs > 0) { int all_scalars; int max_scalar_kind = -1; @@ -1815,7 +1791,7 @@ PyArray_ResultType(npy_intp narrs, PyArrayObject **arr, all_scalars = (ndtypes > 0) ? 0 : 1; /* Compute the maximum "kinds" and whether everything is scalar */ - for (i = 0; i < narrs; ++i) { + for (npy_intp i = 0; i < narrs; ++i) { if (PyArray_NDIM(arr[i]) == 0) { int kind = dtype_kind_to_simplified_ordering( PyArray_DESCR(arr[i])->kind); @@ -1836,7 +1812,7 @@ PyArray_ResultType(npy_intp narrs, PyArrayObject **arr, * If the max scalar kind is bigger than the max array kind, * finish computing the max array kind */ - for (i = 0; i < ndtypes; ++i) { + for (npy_intp i = 0; i < ndtypes; ++i) { int kind = dtype_kind_to_simplified_ordering(dtypes[i]->kind); if (kind > max_array_kind) { max_array_kind = kind; @@ -1848,6 +1824,44 @@ PyArray_ResultType(npy_intp narrs, PyArrayObject **arr, use_min_scalar = 1; } } + return use_min_scalar; +} + + +/*NUMPY_API + * Produces the result type of a bunch of inputs, using the UFunc + * type promotion rules. Use this function when you have a set of + * input arrays, and need to determine an output array dtype. + * + * If all the inputs are scalars (have 0 dimensions) or the maximum "kind" + * of the scalars is greater than the maximum "kind" of the arrays, does + * a regular type promotion. + * + * Otherwise, does a type promotion on the MinScalarType + * of all the inputs. Data types passed directly are treated as array + * types. + * + */ +NPY_NO_EXPORT PyArray_Descr * +PyArray_ResultType(npy_intp narrs, PyArrayObject **arr, + npy_intp ndtypes, PyArray_Descr **dtypes) +{ + npy_intp i; + + /* If there's just one type, pass it through */ + if (narrs + ndtypes == 1) { + PyArray_Descr *ret = NULL; + if (narrs == 1) { + ret = PyArray_DESCR(arr[0]); + } + else { + ret = dtypes[0]; + } + Py_INCREF(ret); + return ret; + } + + int use_min_scalar = should_use_min_scalar(narrs, arr, ndtypes, dtypes); /* Loop through all the types, promoting them */ if (!use_min_scalar) { diff --git a/numpy/core/src/multiarray/convert_datatype.h b/numpy/core/src/multiarray/convert_datatype.h index 653557161..72867ead8 100644 --- a/numpy/core/src/multiarray/convert_datatype.h +++ b/numpy/core/src/multiarray/convert_datatype.h @@ -18,6 +18,10 @@ NPY_NO_EXPORT npy_bool can_cast_scalar_to(PyArray_Descr *scal_type, char *scal_data, PyArray_Descr *to, NPY_CASTING casting); +NPY_NO_EXPORT int +should_use_min_scalar(npy_intp narrs, PyArrayObject **arr, + npy_intp ndtypes, PyArray_Descr **dtypes); + /* * This function calls Py_DECREF on flex_dtype, and replaces it with * a new dtype that has been adapted based on the values in data_dtype diff --git a/numpy/core/src/multiarray/dragon4.c b/numpy/core/src/multiarray/dragon4.c index 8d52672e3..1694596e9 100644 --- a/numpy/core/src/multiarray/dragon4.c +++ b/numpy/core/src/multiarray/dragon4.c @@ -1,31 +1,33 @@ /* * Copyright (c) 2014 Ryan Juckett - * http://www.ryanjuckett.com/ * - * This software is provided 'as-is', without any express or implied - * warranty. In no event will the authors be held liable for any damages - * arising from the use of this software. + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to + * deal in the Software without restriction, including without limitation the + * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or + * sell copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: * - * Permission is granted to anyone to use this software for any purpose, - * including commercial applications, and to alter it and redistribute it - * freely, subject to the following restrictions: + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. * - * 1. The origin of this software must not be misrepresented; you must not - * claim that you wrote the original software. If you use this software - * in a product, an acknowledgment in the product documentation would be - * appreciated but is not required. - * - * 2. Altered source versions must be plainly marked as such, and must not be - * misrepresented as being the original software. - * - * 3. This notice may not be removed or altered from any source - * distribution. + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS + * IN THE SOFTWARE. */ /* * This file contains a modified version of Ryan Juckett's Dragon4 - * implementation, which has been ported from C++ to C and which has + * implementation, obtained from http://www.ryanjuckett.com, + * which has been ported from C++ to C and which has * modifications specific to printing floats in numpy. + * + * Ryan Juckett's original code was under the Zlib license; he gave numpy + * permission to include it under the MIT license instead. */ #include "dragon4.h" diff --git a/numpy/core/src/multiarray/dragon4.h b/numpy/core/src/multiarray/dragon4.h index 2b8b4cef4..3a99bde6c 100644 --- a/numpy/core/src/multiarray/dragon4.h +++ b/numpy/core/src/multiarray/dragon4.h @@ -1,31 +1,33 @@ /* * Copyright (c) 2014 Ryan Juckett - * http://www.ryanjuckett.com/ * - * This software is provided 'as-is', without any express or implied - * warranty. In no event will the authors be held liable for any damages - * arising from the use of this software. + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to + * deal in the Software without restriction, including without limitation the + * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or + * sell copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: * - * Permission is granted to anyone to use this software for any purpose, - * including commercial applications, and to alter it and redistribute it - * freely, subject to the following restrictions: + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. * - * 1. The origin of this software must not be misrepresented; you must not - * claim that you wrote the original software. If you use this software - * in a product, an acknowledgment in the product documentation would be - * appreciated but is not required. - * - * 2. Altered source versions must be plainly marked as such, and must not be - * misrepresented as being the original software. - * - * 3. This notice may not be removed or altered from any source - * distribution. + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS + * IN THE SOFTWARE. */ /* * This file contains a modified version of Ryan Juckett's Dragon4 - * implementation, which has been ported from C++ to C and which has + * implementation, obtained from http://www.ryanjuckett.com, + * which has been ported from C++ to C and which has * modifications specific to printing floats in numpy. + * + * Ryan Juckett's original code was under the Zlib license; he gave numpy + * permission to include it under the MIT license instead. */ #ifndef _NPY_DRAGON4_H_ diff --git a/numpy/core/src/npysort/radixsort.c.src b/numpy/core/src/npysort/radixsort.c.src index c1435bd96..c90b06974 100644 --- a/numpy/core/src/npysort/radixsort.c.src +++ b/numpy/core/src/npysort/radixsort.c.src @@ -50,9 +50,10 @@ nth_byte_@suff@(@type@ key, npy_intp l) { radixsort0_@suff@(@type@ *arr, @type@ *aux, npy_intp num) { npy_intp cnt[sizeof(@type@)][1 << 8] = { { 0 } }; - npy_intp i, l; + npy_intp i; + size_t l; @type@ key0 = KEY_OF(arr[0]); - npy_intp ncols = 0; + size_t ncols = 0; npy_ubyte cols[sizeof(@type@)]; for (i = 0; i < num; i++) { @@ -139,9 +140,10 @@ npy_intp* aradixsort0_@suff@(@type@ *arr, npy_intp *aux, npy_intp *tosort, npy_intp num) { npy_intp cnt[sizeof(@type@)][1 << 8] = { { 0 } }; - npy_intp i, l; + npy_intp i; + size_t l; @type@ key0 = KEY_OF(arr[0]); - npy_intp ncols = 0; + size_t ncols = 0; npy_ubyte cols[sizeof(@type@)]; for (i = 0; i < num; i++) { diff --git a/numpy/core/src/umath/ufunc_type_resolution.c b/numpy/core/src/umath/ufunc_type_resolution.c index 58f915c6e..25dd002ac 100644 --- a/numpy/core/src/umath/ufunc_type_resolution.c +++ b/numpy/core/src/umath/ufunc_type_resolution.c @@ -24,6 +24,7 @@ #include "ufunc_type_resolution.h" #include "ufunc_object.h" #include "common.h" +#include "convert_datatype.h" #include "mem_overlap.h" #if defined(HAVE_CBLAS) @@ -1938,73 +1939,6 @@ type_tuple_userloop_type_resolver(PyUFuncObject *self, return 0; } -/* - * Provides an ordering for the dtype 'kind' character codes, to help - * determine when to use the min_scalar_type function. This groups - * 'kind' into boolean, integer, floating point, and everything else. - */ - -static int -dtype_kind_to_simplified_ordering(char kind) -{ - switch (kind) { - /* Boolean kind */ - case 'b': - return 0; - /* Unsigned int kind */ - case 'u': - /* Signed int kind */ - case 'i': - return 1; - /* Float kind */ - case 'f': - /* Complex kind */ - case 'c': - return 2; - /* Anything else */ - default: - return 3; - } -} - -static int -should_use_min_scalar(PyArrayObject **op, int nop) -{ - int i, use_min_scalar, kind; - int all_scalars = 1, max_scalar_kind = -1, max_array_kind = -1; - - /* - * Determine if there are any scalars, and if so, whether - * the maximum "kind" of the scalars surpasses the maximum - * "kind" of the arrays - */ - use_min_scalar = 0; - if (nop > 1) { - for(i = 0; i < nop; ++i) { - kind = dtype_kind_to_simplified_ordering( - PyArray_DESCR(op[i])->kind); - if (PyArray_NDIM(op[i]) == 0) { - if (kind > max_scalar_kind) { - max_scalar_kind = kind; - } - } - else { - all_scalars = 0; - if (kind > max_array_kind) { - max_array_kind = kind; - } - - } - } - - /* Indicate whether to use the min_scalar_type function */ - if (!all_scalars && max_array_kind >= max_scalar_kind) { - use_min_scalar = 1; - } - } - - return use_min_scalar; -} /* * Does a linear search for the best inner loop of the ufunc. @@ -2030,7 +1964,7 @@ linear_search_type_resolver(PyUFuncObject *self, ufunc_name = ufunc_get_name_cstr(self); - use_min_scalar = should_use_min_scalar(op, nin); + use_min_scalar = should_use_min_scalar(nin, op, 0, NULL); /* If the ufunc has userloops, search for them. */ if (self->userloops) { @@ -2139,7 +2073,7 @@ type_tuple_type_resolver(PyUFuncObject *self, ufunc_name = ufunc_get_name_cstr(self); - use_min_scalar = should_use_min_scalar(op, nin); + use_min_scalar = should_use_min_scalar(nin, op, 0, NULL); /* Fill in specified_types from the tuple or string */ if (PyTuple_Check(type_tup)) { diff --git a/numpy/core/tests/test_dtype.py b/numpy/core/tests/test_dtype.py index 8f33a8daf..4031e7362 100644 --- a/numpy/core/tests/test_dtype.py +++ b/numpy/core/tests/test_dtype.py @@ -11,6 +11,7 @@ from numpy.core._rational_tests import rational from numpy.testing import ( assert_, assert_equal, assert_array_equal, assert_raises, HAS_REFCOUNT) from numpy.compat import pickle +from itertools import permutations def assert_dtype_equal(a, b): assert_equal(a, b) @@ -1178,3 +1179,18 @@ class TestFromCTypes(object): self.check(ctypes.c_uint16.__ctype_be__, np.dtype('>u2')) self.check(ctypes.c_uint8.__ctype_le__, np.dtype('u1')) self.check(ctypes.c_uint8.__ctype_be__, np.dtype('u1')) + + all_types = set(np.typecodes['All']) + all_pairs = permutations(all_types, 2) + + @pytest.mark.parametrize("pair", all_pairs) + def test_pairs(self, pair): + """ + Check that np.dtype('x,y') matches [np.dtype('x'), np.dtype('y')] + Example: np.dtype('d,I') -> dtype([('f0', '<f8'), ('f1', '<u4')]) + """ + # gh-5645: check that np.dtype('i,L') can be used + pair_type = np.dtype('{},{}'.format(*pair)) + expected = np.dtype([('f0', pair[0]), ('f1', pair[1])]) + assert_equal(pair_type, expected) + diff --git a/numpy/distutils/misc_util.py b/numpy/distutils/misc_util.py index 795d348ef..89171eede 100644 --- a/numpy/distutils/misc_util.py +++ b/numpy/distutils/misc_util.py @@ -472,7 +472,7 @@ def is_sequence(seq): return True def is_glob_pattern(s): - return is_string(s) and ('*' in s or '?' is s) + return is_string(s) and ('*' in s or '?' in s) def as_list(seq): if is_sequence(seq): diff --git a/numpy/lib/function_base.py b/numpy/lib/function_base.py index 9690c305f..1fcb6137c 100644 --- a/numpy/lib/function_base.py +++ b/numpy/lib/function_base.py @@ -3722,7 +3722,8 @@ def quantile(a, q, axis=None, out=None, overwrite_input=False, interpolation='linear', keepdims=False): """ Compute the q-th quantile of the data along the specified axis. - ..versionadded:: 1.15.0 + + .. versionadded:: 1.15.0 Parameters ---------- diff --git a/numpy/lib/histograms.py b/numpy/lib/histograms.py index 5bcedc7f4..8474bd5d3 100644 --- a/numpy/lib/histograms.py +++ b/numpy/lib/histograms.py @@ -3,6 +3,7 @@ Histogram-related functions """ from __future__ import division, absolute_import, print_function +import contextlib import functools import operator import warnings @@ -922,7 +923,13 @@ def histogram(a, bins=10, range=None, normed=None, weights=None, def _histogramdd_dispatcher(sample, bins=None, range=None, normed=None, weights=None, density=None): - return (sample, bins, weights) + if hasattr(sample, 'shape'): # same condition as used in histogramdd + yield sample + else: + yield from sample + with contextlib.suppress(TypeError): + yield from bins + yield weights @array_function_dispatch(_histogramdd_dispatcher) diff --git a/numpy/lib/nanfunctions.py b/numpy/lib/nanfunctions.py index f497aca63..9a03d0b39 100644 --- a/numpy/lib/nanfunctions.py +++ b/numpy/lib/nanfunctions.py @@ -1257,6 +1257,7 @@ def nanquantile(a, q, axis=None, out=None, overwrite_input=False, Compute the qth quantile of the data along the specified axis, while ignoring nan values. Returns the qth quantile(s) of the array elements. + .. versionadded:: 1.15.0 Parameters diff --git a/numpy/lib/twodim_base.py b/numpy/lib/twodim_base.py index e165c9b02..0b4e3021a 100644 --- a/numpy/lib/twodim_base.py +++ b/numpy/lib/twodim_base.py @@ -565,7 +565,20 @@ def vander(x, N=None, increasing=False): def _histogram2d_dispatcher(x, y, bins=None, range=None, normed=None, weights=None, density=None): - return (x, y, bins, weights) + yield x + yield y + + # This terrible logic is adapted from the checks in histogram2d + try: + N = len(bins) + except TypeError: + N = 1 + if N != 1 and N != 2: + yield from bins # bins=[x, y] + else: + yield bins + + yield weights @array_function_dispatch(_histogram2d_dispatcher) diff --git a/numpy/linalg/linalg.py b/numpy/linalg/linalg.py index f599f8a58..c90e25686 100644 --- a/numpy/linalg/linalg.py +++ b/numpy/linalg/linalg.py @@ -1509,7 +1509,7 @@ def svd(a, full_matrices=True, compute_uv=True, hermitian=False): enabling a more efficient method for finding singular values. Defaults to False. - ..versionadded:: 1.17.0 + .. versionadded:: 1.17.0 Raises ------ @@ -1913,7 +1913,7 @@ def pinv(a, rcond=1e-15, hermitian=False): enabling a more efficient method for finding singular values. Defaults to False. - ..versionadded:: 1.17.0 + .. versionadded:: 1.17.0 Returns ------- diff --git a/numpy/random/__init__.py b/numpy/random/__init__.py index 21c1adddf..968f3361c 100644 --- a/numpy/random/__init__.py +++ b/numpy/random/__init__.py @@ -98,7 +98,6 @@ set_state Set state of generator. BitGenerator Streams that work with Generator --------------------------------------------- --- MT19937 -DSFMT PCG32 PCG64 Philox @@ -165,7 +164,6 @@ __all__ = [ from . import mtrand from .mtrand import * -from .dsfmt import DSFMT from .generator import Generator from .mt19937 import MT19937 from .pcg32 import PCG32 @@ -176,7 +174,7 @@ from .xoshiro256 import Xoshiro256 from .xoshiro512 import Xoshiro512 from .mtrand import RandomState -__all__ += ['Generator', 'DSFMT', 'MT19937', 'Philox', 'PCG64', 'PCG32', +__all__ += ['Generator', 'MT19937', 'Philox', 'PCG64', 'PCG32', 'ThreeFry', 'Xoshiro256', 'Xoshiro512', 'RandomState'] diff --git a/numpy/random/_pickle.py b/numpy/random/_pickle.py index 307ffe61e..a1f8ab34e 100644 --- a/numpy/random/_pickle.py +++ b/numpy/random/_pickle.py @@ -6,12 +6,10 @@ from .pcg64 import PCG64 from .xoshiro256 import Xoshiro256 from .xoshiro512 import Xoshiro512 -from .dsfmt import DSFMT from .generator import Generator from .mt19937 import MT19937 BitGenerators = {'MT19937': MT19937, - 'DSFMT': DSFMT, 'PCG32': PCG32, 'PCG64': PCG64, 'Philox': Philox, diff --git a/numpy/random/dsfmt.pyx b/numpy/random/dsfmt.pyx deleted file mode 100644 index bca166507..000000000 --- a/numpy/random/dsfmt.pyx +++ /dev/null @@ -1,416 +0,0 @@ -import operator -from cpython.pycapsule cimport PyCapsule_New - -try: - from threading import Lock -except ImportError: - from dummy_threading import Lock - -import numpy as np -cimport numpy as np - -from .common cimport * -from .distributions cimport bitgen_t -from .entropy import random_entropy - -__all__ = ['DSFMT'] - -np.import_array() - -DEF DSFMT_MEXP = 19937 -DEF DSFMT_N = 191 # ((DSFMT_MEXP - 128) / 104 + 1) -DEF DSFMT_N_PLUS_1 = 192 # DSFMT_N + 1 -DEF DSFMT_N64 = DSFMT_N * 2 - -cdef extern from "src/dsfmt/dSFMT.h": - - union W128_T: - uint64_t u[2] - uint32_t u32[4] - double d[2] - - ctypedef W128_T w128_t - - struct DSFMT_T: - w128_t status[DSFMT_N_PLUS_1] - int idx - - ctypedef DSFMT_T dsfmt_t - - struct s_dsfmt_state: - dsfmt_t *state - int has_uint32 - uint32_t uinteger - - double *buffered_uniforms - int buffer_loc - - ctypedef s_dsfmt_state dsfmt_state - - double dsfmt_next_double(dsfmt_state *state) nogil - uint64_t dsfmt_next64(dsfmt_state *state) nogil - uint32_t dsfmt_next32(dsfmt_state *state) nogil - uint64_t dsfmt_next_raw(dsfmt_state *state) nogil - - void dsfmt_init_gen_rand(dsfmt_t *dsfmt, uint32_t seed) - void dsfmt_init_by_array(dsfmt_t *dsfmt, uint32_t init_key[], int key_length) - void dsfmt_jump(dsfmt_state *state) - -cdef uint64_t dsfmt_uint64(void* st) nogil: - return dsfmt_next64(<dsfmt_state *>st) - -cdef uint32_t dsfmt_uint32(void *st) nogil: - return dsfmt_next32(<dsfmt_state *> st) - -cdef double dsfmt_double(void* st) nogil: - return dsfmt_next_double(<dsfmt_state *>st) - -cdef uint64_t dsfmt_raw(void *st) nogil: - return dsfmt_next_raw(<dsfmt_state *>st) - -cdef class DSFMT: - u""" - DSFMT(seed=None) - - Container for the SIMD-based Mersenne Twister pseudo RNG. - - Parameters - ---------- - seed : {None, int, array_like}, optional - Random seed used to initialize the pseudo-random number generator. Can - be any integer between 0 and 2**32 - 1 inclusive, an array (or other - sequence) of unsigned 32-bit integers, or ``None`` (the default). If - `seed` is ``None``, a 32-bit unsigned integer is read from - ``/dev/urandom`` (or the Windows analog) if available. If unavailable, - a 32-bit hash of the time and process ID is used. - - Attributes - ---------- - lock: threading.Lock - Lock instance that is shared so that the same bit git generator can - be used in multiple Generators without corrupting the state. Code that - generates values from a bit generator should hold the bit generator's - lock. - - Notes - ----- - ``DSFMT`` provides a capsule containing function pointers that produce - doubles, and unsigned 32 and 64- bit integers [1]_ . These are not - directly consumable in Python and must be consumed by a ``Generator`` - or similar object that supports low-level access. - - The Python stdlib module "random" also contains a Mersenne Twister - pseudo-random number generator. - - **State and Seeding** - - The ``DSFMT`` state vector consists of a 384 element array of 64-bit - unsigned integers plus a single integer value between 0 and 382 - indicating the current position within the main array. The implementation - used here augments this with a 382 element array of doubles which are used - to efficiently access the random numbers produced by the dSFMT generator. - - ``DSFMT`` is seeded using either a single 32-bit unsigned integer or a - vector of 32-bit unsigned integers. In either case, the input seed is used - as an input (or inputs) for a hashing function, and the output of the - hashing function is used as the initial state. Using a single 32-bit value - for the seed can only initialize a small range of the possible initial - state values. - - **Parallel Features** - - ``DSFMT`` can be used in parallel applications by calling the method - ``jumped`` which advances the state as-if :math:`2^{128}` random numbers - have been generated [2]_. This allows the original sequence to be split - so that distinct segments can be used in each worker process. All - generators should be chained to ensure that the segments come from the same - sequence. - - >>> from numpy.random.entropy import random_entropy - >>> from numpy.random import Generator, DSFMT - >>> seed = random_entropy() - >>> bit_generator = DSFMT(seed) - >>> rg = [] - >>> for _ in range(10): - ... rg.append(Generator(bit_generator)) - ... # Chain the BitGenerators - ... bit_generator = bit_generator.jumped() - - **Compatibility Guarantee** - - ``DSFMT`` makes a guarantee that a fixed seed and will always produce - the same random integer stream. - - References - ---------- - .. [1] Mutsuo Saito and Makoto Matsumoto, "SIMD-oriented Fast Mersenne - Twister: a 128-bit Pseudorandom Number Generator." Monte Carlo - and Quasi-Monte Carlo Methods 2006, Springer, pp. 607--622, 2008. - .. [2] Hiroshi Haramoto, Makoto Matsumoto, and Pierre L\'Ecuyer, "A Fast - Jump Ahead Algorithm for Linear Recurrences in a Polynomial Space", - Sequences and Their Applications - SETA, 290--298, 2008. - """ - cdef dsfmt_state rng_state - cdef bitgen_t _bitgen - cdef public object capsule - cdef object _cffi - cdef object _ctypes - cdef public object lock - - def __init__(self, seed=None): - self.rng_state.state = <dsfmt_t *>PyArray_malloc_aligned(sizeof(dsfmt_t)) - self.rng_state.buffered_uniforms = <double *>PyArray_calloc_aligned(DSFMT_N64, sizeof(double)) - self.rng_state.buffer_loc = DSFMT_N64 - self.seed(seed) - self.lock = Lock() - - self._bitgen.state = <void *>&self.rng_state - self._bitgen.next_uint64 = &dsfmt_uint64 - self._bitgen.next_uint32 = &dsfmt_uint32 - self._bitgen.next_double = &dsfmt_double - self._bitgen.next_raw = &dsfmt_raw - cdef const char *name = "BitGenerator" - self.capsule = PyCapsule_New(<void *>&self._bitgen, name, NULL) - - self._cffi = None - self._ctypes = None - - # Pickling support: - def __getstate__(self): - return self.state - - def __setstate__(self, state): - self.state = state - - def __reduce__(self): - from ._pickle import __bit_generator_ctor - return __bit_generator_ctor, (self.state['bit_generator'],), self.state - - def __dealloc__(self): - if self.rng_state.state: - PyArray_free_aligned(self.rng_state.state) - if self.rng_state.buffered_uniforms: - PyArray_free_aligned(self.rng_state.buffered_uniforms) - - cdef _reset_state_variables(self): - self.rng_state.buffer_loc = DSFMT_N64 - - def random_raw(self, size=None, output=True): - """ - random_raw(self, size=None) - - Return randoms as generated by the underlying BitGenerator - - Parameters - ---------- - size : int or tuple of ints, optional - 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. - output : bool, optional - Output values. Used for performance testing since the generated - values are not returned. - - Returns - ------- - out : uint or ndarray - Drawn samples. - - Notes - ----- - This method directly exposes the the raw underlying pseudo-random - number generator. All values are returned as unsigned 64-bit - values irrespective of the number of bits produced by the PRNG. - - See the class docstring for the number of bits returned. - """ - return random_raw(&self._bitgen, self.lock, size, output) - - def _benchmark(self, Py_ssize_t cnt, method=u'uint64'): - return benchmark(&self._bitgen, self.lock, cnt, method) - - def seed(self, seed=None): - """ - seed(seed=None) - - Seed the generator. - - Parameters - ---------- - seed : {None, int, array_like}, optional - Random seed initializing the pseudo-random number generator. - Can be an integer in [0, 2**32-1], array of integers in - [0, 2**32-1] or ``None`` (the default). If `seed` is ``None``, - then ``DSFMT`` will try to read entropy from ``/dev/urandom`` - (or the Windows analog) if available to produce a 32-bit - seed. If unavailable, a 32-bit hash of the time and process - ID is used. - - Raises - ------ - ValueError - If seed values are out of range for the PRNG. - """ - cdef np.ndarray obj - try: - if seed is None: - try: - seed = random_entropy(1) - except RuntimeError: - seed = random_entropy(1, 'fallback') - dsfmt_init_gen_rand(self.rng_state.state, seed) - else: - if hasattr(seed, 'squeeze'): - seed = seed.squeeze() - idx = operator.index(seed) - if idx > int(2**32 - 1) or idx < 0: - raise ValueError("Seed must be between 0 and 2**32 - 1") - dsfmt_init_gen_rand(self.rng_state.state, seed) - except TypeError: - obj = np.asarray(seed).astype(np.int64, casting='safe').ravel() - if ((obj > int(2**32 - 1)) | (obj < 0)).any(): - raise ValueError("Seed must be between 0 and 2**32 - 1") - obj = obj.astype(np.uint32, casting='unsafe', order='C') - dsfmt_init_by_array(self.rng_state.state, - <uint32_t *>obj.data, - np.PyArray_DIM(obj, 0)) - # Clear the buffer - self._reset_state_variables() - - cdef jump_inplace(self, iter): - """ - Jump state in-place - - Not part of public API - - Parameters - ---------- - iter : integer, positive - Number of times to jump the state of the rng. - """ - cdef np.npy_intp i - for i in range(iter): - dsfmt_jump(&self.rng_state) - # Clear the buffer - self._reset_state_variables() - - def jumped(self, np.npy_intp jumps=1): - """ - jumped(jumps=1) - - Returns a new bit generator with the state jumped - - The state of the returned big generator is jumped as-if - 2**(128 * jumps) random numbers have been generated. - - Parameters - ---------- - jumps : integer, positive - Number of times to jump the state of the bit generator returned - - Returns - ------- - bit_generator : DSFMT - New instance of generator jumped iter times - """ - cdef DSFMT bit_generator - - bit_generator = self.__class__() - bit_generator.state = self.state - bit_generator.jump_inplace(jumps) - - return bit_generator - - @property - def state(self): - """ - Get or set the PRNG state - - Returns - ------- - state : dict - Dictionary containing the information required to describe the - state of the PRNG - """ - - cdef Py_ssize_t i, j, loc = 0 - cdef uint64_t[::1] state - cdef double[::1] buffered_uniforms - - state = np.empty(2 *DSFMT_N_PLUS_1, dtype=np.uint64) - for i in range(DSFMT_N_PLUS_1): - for j in range(2): - state[loc] = self.rng_state.state.status[i].u[j] - loc += 1 - buffered_uniforms = np.empty(DSFMT_N64, dtype=np.double) - for i in range(DSFMT_N64): - buffered_uniforms[i] = self.rng_state.buffered_uniforms[i] - return {'bit_generator': self.__class__.__name__, - 'state': {'state': np.asarray(state), - 'idx': self.rng_state.state.idx}, - 'buffer_loc': self.rng_state.buffer_loc, - 'buffered_uniforms': np.asarray(buffered_uniforms)} - - @state.setter - def state(self, value): - cdef Py_ssize_t i, j, loc = 0 - if not isinstance(value, dict): - raise TypeError('state must be a dict') - bitgen = value.get('bit_generator', '') - if bitgen != self.__class__.__name__: - raise ValueError('state must be for a {0} ' - 'PRNG'.format(self.__class__.__name__)) - state = value['state']['state'] - for i in range(DSFMT_N_PLUS_1): - for j in range(2): - self.rng_state.state.status[i].u[j] = state[loc] - loc += 1 - self.rng_state.state.idx = value['state']['idx'] - buffered_uniforms = value['buffered_uniforms'] - for i in range(DSFMT_N64): - self.rng_state.buffered_uniforms[i] = buffered_uniforms[i] - self.rng_state.buffer_loc = value['buffer_loc'] - - @property - def ctypes(self): - """ - ctypes interface - - Returns - ------- - interface : namedtuple - Named tuple containing ctypes wrapper - - * state_address - Memory address of the state struct - * state - pointer to the state struct - * next_uint64 - function pointer to produce 64 bit integers - * next_uint32 - function pointer to produce 32 bit integers - * next_double - function pointer to produce doubles - * bitgen - pointer to the bit generator struct - """ - if self._ctypes is None: - self._ctypes = prepare_ctypes(&self._bitgen) - - return self._ctypes - - @property - def cffi(self): - """ - CFFI interface - - Returns - ------- - interface : namedtuple - Named tuple containing CFFI wrapper - - * state_address - Memory address of the state struct - * state - pointer to the state struct - * next_uint64 - function pointer to produce 64 bit integers - * next_uint32 - function pointer to produce 32 bit integers - * next_double - function pointer to produce doubles - * bitgen - pointer to the bit generator struct - """ - if self._cffi is not None: - return self._cffi - self._cffi = prepare_cffi(&self._bitgen) - return self._cffi diff --git a/numpy/random/generator.pyx b/numpy/random/generator.pyx index 4e2b18f44..2d2979441 100644 --- a/numpy/random/generator.pyx +++ b/numpy/random/generator.pyx @@ -3095,9 +3095,11 @@ cdef class Generator: Parameters ---------- ngood : int or array_like of ints - Number of ways to make a good selection. Must be nonnegative. + Number of ways to make a good selection. Must be nonnegative and + less than 10**9. nbad : int or array_like of ints - Number of ways to make a bad selection. Must be nonnegative. + Number of ways to make a bad selection. Must be nonnegative and + less than 10**9. nsample : int or array_like of ints Number of items sampled. Must be nonnegative and less than ``ngood + nbad``. @@ -3142,6 +3144,13 @@ cdef class Generator: replacement (or the sample space is infinite). As the sample space becomes large, this distribution approaches the binomial. + The arguments `ngood` and `nbad` each must be less than `10**9`. For + extremely large arguments, the algorithm that is used to compute the + samples [4]_ breaks down because of loss of precision in floating point + calculations. For such large values, if `nsample` is not also large, + the distribution can be approximated with the binomial distribution, + `binomial(n=nsample, p=ngood/(ngood + nbad))`. + References ---------- .. [1] Lentner, Marvin, "Elementary Applied Statistics", Bogden @@ -3151,6 +3160,9 @@ cdef class Generator: http://mathworld.wolfram.com/HypergeometricDistribution.html .. [3] Wikipedia, "Hypergeometric distribution", https://en.wikipedia.org/wiki/Hypergeometric_distribution + .. [4] Stadlober, Ernst, "The ratio of uniforms approach for generating + discrete random variates", Journal of Computational and Applied + Mathematics, 31, pp. 181-189 (1990). Examples -------- @@ -3172,6 +3184,7 @@ cdef class Generator: # answer = 0.003 ... pretty unlikely! """ + DEF HYPERGEOM_MAX = 10**9 cdef bint is_scalar = True cdef np.ndarray ongood, onbad, onsample cdef int64_t lngood, lnbad, lnsample @@ -3186,6 +3199,9 @@ cdef class Generator: lnbad = <int64_t>nbad lnsample = <int64_t>nsample + if lngood >= HYPERGEOM_MAX or lnbad >= HYPERGEOM_MAX: + raise ValueError("both ngood and nbad must be less than %d" % + HYPERGEOM_MAX) if lngood + lnbad < lnsample: raise ValueError("ngood + nbad < nsample") return disc(&random_hypergeometric, &self._bitgen, size, self.lock, 0, 3, @@ -3193,8 +3209,13 @@ cdef class Generator: lnbad, 'nbad', CONS_NON_NEGATIVE, lnsample, 'nsample', CONS_NON_NEGATIVE) + if np.any(ongood >= HYPERGEOM_MAX) or np.any(onbad >= HYPERGEOM_MAX): + raise ValueError("both ngood and nbad must be less than %d" % + HYPERGEOM_MAX) + if np.any(np.less(np.add(ongood, onbad), onsample)): raise ValueError("ngood + nbad < nsample") + return discrete_broadcast_iii(&random_hypergeometric, &self._bitgen, size, self.lock, ongood, 'ngood', CONS_NON_NEGATIVE, onbad, 'nbad', CONS_NON_NEGATIVE, diff --git a/numpy/random/setup.py b/numpy/random/setup.py index 5f467a2f5..84e9f1c68 100644 --- a/numpy/random/setup.py +++ b/numpy/random/setup.py @@ -55,9 +55,6 @@ def configuration(parent_package='', top_path=None): # Use legacy integer variable sizes LEGACY_DEFS = [('NP_RANDOM_LEGACY', '1')] - # Required defined for DSFMT size and to allow it to detect SSE2 using - # config file information - DSFMT_DEFS = [('DSFMT_MEXP', '19937'), ("HAVE_NPY_CONFIG_H", "1")] PCG64_DEFS = [] if 1 or sys.maxsize < 2 ** 32 or os.name == 'nt': # Force emulated mode here @@ -75,19 +72,6 @@ def configuration(parent_package='', top_path=None): ], define_macros=defs, ) - config.add_extension('dsfmt', - sources=['dsfmt.c', 'src/dsfmt/dSFMT.c', - 'src/dsfmt/dSFMT-jump.c', - 'src/aligned_malloc/aligned_malloc.c'], - include_dirs=['.', 'src', join('src', 'dsfmt')], - libraries=EXTRA_LIBRARIES, - extra_compile_args=EXTRA_COMPILE_ARGS, - extra_link_args=EXTRA_LINK_ARGS, - depends=[join('src', 'dsfmt', 'dsfmt.h'), - 'dsfmt.pyx', - ], - define_macros=defs + DSFMT_DEFS, - ) for gen in ['mt19937']: # gen.pyx, src/gen/gen.c, src/gen/gen-jump.c config.add_extension(gen, @@ -126,12 +110,15 @@ def configuration(parent_package='', top_path=None): depends=['%s.pyx' % gen], define_macros=defs, ) + other_srcs = [ + 'src/distributions/logfactorial.c', + 'src/distributions/distributions.c', + 'src/distributions/random_hypergeometric.c', + ] for gen in ['generator', 'bounded_integers']: # gen.pyx, src/distributions/distributions.c config.add_extension(gen, - sources=['{0}.c'.format(gen), - join('src', 'distributions', - 'distributions.c')], + sources=['{0}.c'.format(gen)] + other_srcs, libraries=EXTRA_LIBRARIES, extra_compile_args=EXTRA_COMPILE_ARGS, include_dirs=['.', 'src'], @@ -140,15 +127,17 @@ def configuration(parent_package='', top_path=None): define_macros=defs, ) config.add_extension('mtrand', + # mtrand does not depend on random_hypergeometric.c. sources=['mtrand.c', 'src/legacy/legacy-distributions.c', + 'src/distributions/logfactorial.c', 'src/distributions/distributions.c'], include_dirs=['.', 'src', 'src/legacy'], libraries=EXTRA_LIBRARIES, extra_compile_args=EXTRA_COMPILE_ARGS, extra_link_args=EXTRA_LINK_ARGS, depends=['mtrand.pyx'], - define_macros=defs + DSFMT_DEFS + LEGACY_DEFS, + define_macros=defs + LEGACY_DEFS, ) return config diff --git a/numpy/random/src/distributions/distributions.c b/numpy/random/src/distributions/distributions.c index def29d850..65257ecbf 100644 --- a/numpy/random/src/distributions/distributions.c +++ b/numpy/random/src/distributions/distributions.c @@ -1,5 +1,6 @@ #include "distributions.h" #include "ziggurat_constants.h" +#include "logfactorial.h" #if defined(_MSC_VER) && defined(_WIN64) #include <intrin.h> @@ -169,7 +170,7 @@ float random_standard_exponential_zig_f(bitgen_t *bitgen_state) { static NPY_INLINE double next_gauss_zig(bitgen_t *bitgen_state) { uint64_t r; int sign; - int64_t rabs; + uint64_t rabs; int idx; double x, xx, yy; for (;;) { @@ -178,7 +179,7 @@ static NPY_INLINE double next_gauss_zig(bitgen_t *bitgen_state) { idx = r & 0xff; r >>= 8; sign = r & 0x1; - rabs = (int64_t)((r >> 1) & 0x000fffffffffffff); + rabs = (r >> 1) & 0x000fffffffffffff; x = rabs * wi_double[idx]; if (sign & 0x1) x = -x; @@ -215,7 +216,7 @@ void random_gauss_zig_fill(bitgen_t *bitgen_state, npy_intp cnt, double *out) { float random_gauss_zig_f(bitgen_t *bitgen_state) { uint32_t r; int sign; - int32_t rabs; + uint32_t rabs; int idx; float x, xx, yy; for (;;) { @@ -223,7 +224,7 @@ float random_gauss_zig_f(bitgen_t *bitgen_state) { r = next_uint32(bitgen_state); idx = r & 0xff; sign = (r >> 8) & 0x1; - rabs = (int32_t)((r >> 9) & 0x0007fffff); + rabs = (r >> 9) & 0x0007fffff; x = rabs * wi_float[idx]; if (sign & 0x1) x = -x; @@ -468,8 +469,11 @@ uint64_t random_uint(bitgen_t *bitgen_state) { * log-gamma function to support some of these distributions. The * algorithm comes from SPECFUN by Shanjie Zhang and Jianming Jin and their * book "Computation of Special Functions", 1996, John Wiley & Sons, Inc. + * + * If loggam(k+1) is being used to compute log(k!) for an integer k, consider + * using logfactorial(k) instead. */ -static double loggam(double x) { +double loggam(double x) { double x0, x2, xp, gl, gl0; RAND_INT_TYPE k, n; @@ -1127,105 +1131,6 @@ double random_triangular(bitgen_t *bitgen_state, double left, double mode, } } -RAND_INT_TYPE random_hypergeometric_hyp(bitgen_t *bitgen_state, - RAND_INT_TYPE good, RAND_INT_TYPE bad, - RAND_INT_TYPE sample) { - RAND_INT_TYPE d1, k, z; - double d2, u, y; - - d1 = bad + good - sample; - d2 = (double)MIN(bad, good); - - y = d2; - k = sample; - while (y > 0.0) { - u = next_double(bitgen_state); - y -= (RAND_INT_TYPE)floor(u + y / (d1 + k)); - k--; - if (k == 0) - break; - } - z = (RAND_INT_TYPE)(d2 - y); - if (good > bad) - z = sample - z; - return z; -} - -/* D1 = 2*sqrt(2/e) */ -/* D2 = 3 - 2*sqrt(3/e) */ -#define D1 1.7155277699214135 -#define D2 0.8989161620588988 -RAND_INT_TYPE random_hypergeometric_hrua(bitgen_t *bitgen_state, - RAND_INT_TYPE good, RAND_INT_TYPE bad, - RAND_INT_TYPE sample) { - RAND_INT_TYPE mingoodbad, maxgoodbad, popsize, m, d9; - double d4, d5, d6, d7, d8, d10, d11; - RAND_INT_TYPE Z; - double T, W, X, Y; - - mingoodbad = MIN(good, bad); - popsize = good + bad; - maxgoodbad = MAX(good, bad); - m = MIN(sample, popsize - sample); - d4 = ((double)mingoodbad) / popsize; - d5 = 1.0 - d4; - d6 = m * d4 + 0.5; - d7 = sqrt((double)(popsize - m) * sample * d4 * d5 / (popsize - 1) + 0.5); - d8 = D1 * d7 + D2; - d9 = (RAND_INT_TYPE)floor((double)(m + 1) * (mingoodbad + 1) / (popsize + 2)); - d10 = (loggam(d9 + 1) + loggam(mingoodbad - d9 + 1) + loggam(m - d9 + 1) + - loggam(maxgoodbad - m + d9 + 1)); - d11 = MIN(MIN(m, mingoodbad) + 1.0, floor(d6 + 16 * d7)); - /* 16 for 16-decimal-digit precision in D1 and D2 */ - - while (1) { - X = next_double(bitgen_state); - Y = next_double(bitgen_state); - W = d6 + d8 * (Y - 0.5) / X; - - /* fast rejection: */ - if ((W < 0.0) || (W >= d11)) - continue; - - Z = (RAND_INT_TYPE)floor(W); - T = d10 - (loggam(Z + 1) + loggam(mingoodbad - Z + 1) + loggam(m - Z + 1) + - loggam(maxgoodbad - m + Z + 1)); - - /* fast acceptance: */ - if ((X * (4.0 - X) - 3.0) <= T) - break; - - /* fast rejection: */ - if (X * (X - T) >= 1) - continue; - /* log(0.0) is ok here, since always accept */ - if (2.0 * log(X) <= T) - break; /* acceptance */ - } - - /* this is a correction to HRUA* by Ivan Frohne in rv.py */ - if (good > bad) - Z = m - Z; - - /* another fix from rv.py to allow sample to exceed popsize/2 */ - if (m < sample) - Z = good - Z; - - return Z; -} -#undef D1 -#undef D2 - -RAND_INT_TYPE random_hypergeometric(bitgen_t *bitgen_state, RAND_INT_TYPE good, - RAND_INT_TYPE bad, RAND_INT_TYPE sample) { - if (sample > 10) { - return random_hypergeometric_hrua(bitgen_state, good, bad, sample); - } else if (sample > 0) { - return random_hypergeometric_hyp(bitgen_state, good, bad, sample); - } else { - return 0; - } -} uint64_t random_interval(bitgen_t *bitgen_state, uint64_t max) { uint64_t mask, value; diff --git a/numpy/random/src/distributions/distributions.h b/numpy/random/src/distributions/distributions.h index e43a5279c..3178725af 100644 --- a/numpy/random/src/distributions/distributions.h +++ b/numpy/random/src/distributions/distributions.h @@ -84,6 +84,8 @@ static NPY_INLINE double next_double(bitgen_t *bitgen_state) { return bitgen_state->next_double(bitgen_state->state); } +DECLDIR double loggam(double x); + DECLDIR float random_float(bitgen_t *bitgen_state); DECLDIR double random_double(bitgen_t *bitgen_state); DECLDIR void random_double_fill(bitgen_t *bitgen_state, npy_intp cnt, double *out); @@ -160,8 +162,8 @@ DECLDIR RAND_INT_TYPE random_geometric_search(bitgen_t *bitgen_state, double p); DECLDIR RAND_INT_TYPE random_geometric_inversion(bitgen_t *bitgen_state, double p); DECLDIR RAND_INT_TYPE random_geometric(bitgen_t *bitgen_state, double p); DECLDIR RAND_INT_TYPE random_zipf(bitgen_t *bitgen_state, double a); -DECLDIR RAND_INT_TYPE random_hypergeometric(bitgen_t *bitgen_state, RAND_INT_TYPE good, - RAND_INT_TYPE bad, RAND_INT_TYPE sample); +DECLDIR int64_t random_hypergeometric(bitgen_t *bitgen_state, + int64_t good, int64_t bad, int64_t sample); DECLDIR uint64_t random_interval(bitgen_t *bitgen_state, uint64_t max); diff --git a/numpy/random/src/distributions/logfactorial.c b/numpy/random/src/distributions/logfactorial.c new file mode 100644 index 000000000..130516469 --- /dev/null +++ b/numpy/random/src/distributions/logfactorial.c @@ -0,0 +1,158 @@ + +#include <math.h> +#include <stdint.h> + +/* + * logfact[k] holds log(k!) for k = 0, 1, 2, ..., 125. + */ + +static const double logfact[] = { + 0, + 0, + 0.69314718055994529, + 1.791759469228055, + 3.1780538303479458, + 4.7874917427820458, + 6.5792512120101012, + 8.5251613610654147, + 10.604602902745251, + 12.801827480081469, + 15.104412573075516, + 17.502307845873887, + 19.987214495661885, + 22.552163853123425, + 25.19122118273868, + 27.89927138384089, + 30.671860106080672, + 33.505073450136891, + 36.395445208033053, + 39.339884187199495, + 42.335616460753485, + 45.380138898476908, + 48.471181351835227, + 51.606675567764377, + 54.784729398112319, + 58.003605222980518, + 61.261701761002001, + 64.557538627006338, + 67.88974313718154, + 71.257038967168015, + 74.658236348830158, + 78.092223553315307, + 81.557959456115043, + 85.054467017581516, + 88.580827542197682, + 92.136175603687093, + 95.719694542143202, + 99.330612454787428, + 102.96819861451381, + 106.63176026064346, + 110.32063971475739, + 114.03421178146171, + 117.77188139974507, + 121.53308151543864, + 125.3172711493569, + 129.12393363912722, + 132.95257503561632, + 136.80272263732635, + 140.67392364823425, + 144.5657439463449, + 148.47776695177302, + 152.40959258449735, + 156.3608363030788, + 160.3311282166309, + 164.32011226319517, + 168.32744544842765, + 172.35279713916279, + 176.39584840699735, + 180.45629141754378, + 184.53382886144948, + 188.6281734236716, + 192.7390472878449, + 196.86618167289001, + 201.00931639928152, + 205.1681994826412, + 209.34258675253685, + 213.53224149456327, + 217.73693411395422, + 221.95644181913033, + 226.1905483237276, + 230.43904356577696, + 234.70172344281826, + 238.97838956183432, + 243.26884900298271, + 247.57291409618688, + 251.89040220972319, + 256.22113555000954, + 260.56494097186322, + 264.92164979855278, + 269.29109765101981, + 273.67312428569369, + 278.06757344036612, + 282.4742926876304, + 286.89313329542699, + 291.32395009427029, + 295.76660135076065, + 300.22094864701415, + 304.68685676566872, + 309.1641935801469, + 313.65282994987905, + 318.1526396202093, + 322.66349912672615, + 327.1852877037752, + 331.71788719692847, + 336.26118197919845, + 340.81505887079902, + 345.37940706226686, + 349.95411804077025, + 354.53908551944079, + 359.1342053695754, + 363.73937555556347, + 368.35449607240474, + 372.97946888568902, + 377.61419787391867, + 382.25858877306001, + 386.91254912321756, + 391.57598821732961, + 396.24881705179155, + 400.93094827891576, + 405.6222961611449, + 410.32277652693733, + 415.03230672824964, + 419.75080559954472, + 424.47819341825709, + 429.21439186665157, + 433.95932399501481, + 438.71291418612117, + 443.47508812091894, + 448.24577274538461, + 453.02489623849613, + 457.81238798127816, + 462.60817852687489, + 467.4121995716082, + 472.22438392698058, + 477.04466549258564, + 481.87297922988796 +}; + +/* + * Compute log(k!) + */ + +double logfactorial(int64_t k) +{ + const double halfln2pi = 0.9189385332046728; + + if (k < (int64_t) (sizeof(logfact)/sizeof(logfact[0]))) { + /* Use the lookup table. */ + return logfact[k]; + } + + /* + * Use the Stirling series, truncated at the 1/k**3 term. + * (In a Python implementation of this approximation, the result + * was within 2 ULP of the best 64 bit floating point value for + * k up to 10000000.) + */ + return (k + 0.5)*log(k) - k + (halfln2pi + (1.0/k)*(1/12.0 - 1/(360.0*k*k))); +} diff --git a/numpy/random/src/distributions/logfactorial.h b/numpy/random/src/distributions/logfactorial.h new file mode 100644 index 000000000..1fedef3f6 --- /dev/null +++ b/numpy/random/src/distributions/logfactorial.h @@ -0,0 +1,9 @@ + +#ifndef LOGFACTORIAL_H +#define LOGFACTORIAL_H + +#include <stdint.h> + +double logfactorial(int64_t k); + +#endif diff --git a/numpy/random/src/distributions/random_hypergeometric.c b/numpy/random/src/distributions/random_hypergeometric.c new file mode 100644 index 000000000..59a3a4b9b --- /dev/null +++ b/numpy/random/src/distributions/random_hypergeometric.c @@ -0,0 +1,260 @@ +#include <stdint.h> +#include "distributions.h" +#include "logfactorial.h" + +/* + * Generate a sample from the hypergeometric distribution. + * + * Assume sample is not greater than half the total. See below + * for how the opposite case is handled. + * + * We initialize the following: + * computed_sample = sample + * remaining_good = good + * remaining_total = good + bad + * + * In the loop: + * * computed_sample counts down to 0; + * * remaining_good is the number of good choices not selected yet; + * * remaining_total is the total number of choices not selected yet. + * + * In the loop, we select items by choosing a random integer in + * the interval [0, remaining_total), and if the value is less + * than remaining_good, it means we have selected a good one, + * so remaining_good is decremented. Then, regardless of that + * result, computed_sample is decremented. The loop continues + * until either computed_sample is 0, remaining_good is 0, or + * remaining_total == remaining_good. In the latter case, it + * means there are only good choices left, so we can stop the + * loop early and select what is left of computed_sample from + * the good choices (i.e. decrease remaining_good by computed_sample). + * + * When the loop exits, the actual number of good choices is + * good - remaining_good. + * + * If sample is more than half the total, then initially we set + * computed_sample = total - sample + * and at the end we return remaining_good (i.e. the loop in effect + * selects the complement of the result). + * + * It is assumed that when this function is called: + * * good, bad and sample are nonnegative; + * * the sum good+bad will not result in overflow; + * * sample <= good+bad. + */ + +static int64_t hypergeometric_sample(bitgen_t *bitgen_state, + int64_t good, int64_t bad, int64_t sample) +{ + int64_t remaining_total, remaining_good, result, computed_sample; + int64_t total = good + bad; + + if (sample > total/2) { + computed_sample = total - sample; + } + else { + computed_sample = sample; + } + + remaining_total = total; + remaining_good = good; + + while ((computed_sample > 0) && (remaining_good > 0) && + (remaining_total > remaining_good)) { + // random_interval(bitgen_state, max) returns an integer in + // [0, max] *inclusive*, so we decrement remaining_total before + // passing it to random_interval(). + --remaining_total; + if ((int64_t) random_interval(bitgen_state, + remaining_total) < remaining_good) { + // Selected a "good" one, so decrement remaining_good. + --remaining_good; + } + --computed_sample; + } + + if (remaining_total == remaining_good) { + // Only "good" choices are left. + remaining_good -= computed_sample; + } + + if (sample > total/2) { + result = remaining_good; + } + else { + result = good - remaining_good; + } + + return result; +} + + +// D1 = 2*sqrt(2/e) +// D2 = 3 - 2*sqrt(3/e) +#define D1 1.7155277699214135 +#define D2 0.8989161620588988 + +/* + * Generate variates from the hypergeometric distribution + * using the ratio-of-uniforms method. + * + * In the code, the variable names a, b, c, g, h, m, p, q, K, T, + * U and X match the names used in "Algorithm HRUA" beginning on + * page 82 of Stadlober's 1989 thesis. + * + * It is assumed that when this function is called: + * * good, bad and sample are nonnegative; + * * the sum good+bad will not result in overflow; + * * sample <= good+bad. + * + * References: + * - Ernst Stadlober's thesis "Sampling from Poisson, Binomial and + * Hypergeometric Distributions: Ratio of Uniforms as a Simple and + * Fast Alternative" (1989) + * - Ernst Stadlober, "The ratio of uniforms approach for generating + * discrete random variates", Journal of Computational and Applied + * Mathematics, 31, pp. 181-189 (1990). + */ + +static int64_t hypergeometric_hrua(bitgen_t *bitgen_state, + int64_t good, int64_t bad, int64_t sample) +{ + int64_t mingoodbad, maxgoodbad, popsize; + int64_t computed_sample; + double p, q; + double mu, var; + double a, c, b, h, g; + int64_t m, K; + + popsize = good + bad; + computed_sample = MIN(sample, popsize - sample); + mingoodbad = MIN(good, bad); + maxgoodbad = MAX(good, bad); + + /* + * Variables that do not match Stadlober (1989) + * Here Stadlober + * ---------------- --------- + * mingoodbad M + * popsize N + * computed_sample n + */ + + p = ((double) mingoodbad) / popsize; + q = ((double) maxgoodbad) / popsize; + + // mu is the mean of the distribution. + mu = computed_sample * p; + + a = mu + 0.5; + + // var is the variance of the distribution. + var = ((double)(popsize - computed_sample) * + computed_sample * p * q / (popsize - 1)); + + c = sqrt(var + 0.5); + + /* + * h is 2*s_hat (See Stadlober's theses (1989), Eq. (5.17); or + * Stadlober (1990), Eq. 8). s_hat is the scale of the "table mountain" + * function that dominates the scaled hypergeometric PMF ("scaled" means + * normalized to have a maximum value of 1). + */ + h = D1*c + D2; + + m = (int64_t) floor((double)(computed_sample + 1) * (mingoodbad + 1) / + (popsize + 2)); + + g = (logfactorial(m) + + logfactorial(mingoodbad - m) + + logfactorial(computed_sample - m) + + logfactorial(maxgoodbad - computed_sample + m)); + + /* + * b is the upper bound for random samples: + * ... min(computed_sample, mingoodbad) + 1 is the length of the support. + * ... floor(a + 16*c) is 16 standard deviations beyond the mean. + * + * The idea behind the second upper bound is that values that far out in + * the tail have negligible probabilities. + * + * There is a comment in a previous version of this algorithm that says + * "16 for 16-decimal-digit precision in D1 and D2", + * but there is no documented justification for this value. A lower value + * might work just as well, but I've kept the value 16 here. + */ + b = MIN(MIN(computed_sample, mingoodbad) + 1, floor(a + 16*c)); + + while (1) { + double U, V, X, T; + double gp; + U = random_double(bitgen_state); + V = random_double(bitgen_state); // "U star" in Stadlober (1989) + X = a + h*(V - 0.5) / U; + + // fast rejection: + if ((X < 0.0) || (X >= b)) { + continue; + } + + K = (int64_t) floor(X); + + gp = (logfactorial(K) + + logfactorial(mingoodbad - K) + + logfactorial(computed_sample - K) + + logfactorial(maxgoodbad - computed_sample + K)); + + T = g - gp; + + // fast acceptance: + if ((U*(4.0 - U) - 3.0) <= T) { + break; + } + + // fast rejection: + if (U*(U - T) >= 1) { + continue; + } + + if (2.0*log(U) <= T) { + // acceptance + break; + } + } + + if (good > bad) { + K = computed_sample - K; + } + + if (computed_sample < sample) { + K = good - K; + } + + return K; +} + + +/* + * Draw a sample from the hypergeometric distribution. + * + * It is assumed that when this function is called: + * * good, bad and sample are nonnegative; + * * the sum good+bad will not result in overflow; + * * sample <= good+bad. + */ + +int64_t random_hypergeometric(bitgen_t *bitgen_state, + int64_t good, int64_t bad, int64_t sample) +{ + int64_t r; + + if ((sample >= 10) && (sample <= good + bad - 10)) { + // This will use the ratio-of-uniforms method. + r = hypergeometric_hrua(bitgen_state, good, bad, sample); + } + else { + // The simpler implementation is faster for small samples. + r = hypergeometric_sample(bitgen_state, good, bad, sample); + } + return r; +} diff --git a/numpy/random/src/dsfmt/128-bit-jump.poly.txt b/numpy/random/src/dsfmt/128-bit-jump.poly.txt deleted file mode 100644 index fea1318fb..000000000 --- a/numpy/random/src/dsfmt/128-bit-jump.poly.txt +++ /dev/null @@ -1,2 +0,0 @@ -jump polynomial: -f4dfa6c62049d0776e0bf6f1e953f3aa38abb113df86be024eab3773ad5f2b82ead936022e656dff7e562691c59dd5f7d2566b78d9669002503c4ddb1888a49f32333f515e6c60c4ecd221078ec6f26f0a90f4875067ca1f399a99775037adf905566e2c7e6b42131420f8f04f112c92621c9b1502f2a8aefad6c667904af62f0d55e02d396902d3b89450103c5ce5fe0408d97cbb864861b49e4e42048ff3310b48faac55095a7f422eea4aade752f947f947c6be0a0c665bdea099246ab9eff658ea8ca468bf49d0227748367878de06d7bd86ea6708fcac6e252f5f00f04309b2aac3036b64afb39d990427c6c9f03477cc7e935c43c0e61bc161db8eb15516eee8cb377ecbc1849207990fb6778721b29bfe0d89bfda1b3772fa5b0b1f7ec3daf36052032285898c6f6396f55010c31f8201b7e2e51d94f920bfe57684c5415cc342cb39a0045d9793d13cf8646096daeb8bb9bfc20a90de8f2426da8733267a9b9674f32154e8f84a9932223a2ca3c787d0b66df6675febbdfcba2f9cef09c621c57e11098b3289c77397aaae8b104642ffe0c4b75598efbc53745984d68b4d6656cae299ae2be55217a9a02b009ca7be32f47fbe434bce4914a34d0c9b0085bede9b8a99319c34660d66f0124b5a7714c4bf3cbfec3ee43ed817087168bad80133bebaeeb68cf7929a24d1bb3de831a8340d220906ab04159cf94b21d5ee813bd7c80f10f01b43052af530917513b169254c25d6fcfe6cb420d6ce92f54886ef6eaf9a5ba35e893ff593834d05ddf28899e42d729c7df3d21ef036020789739366f0c11ec52ff92a0bfd8ba69508e27b20fabb8217bd36b90e5aa918159ac87913bc7b46c04e366c23c92807fbe9c6a407e6a4db0b4fc23c3b6c706b5ca058fe8c190f849f18d16d6b48b5ed760eb202fd566291a799420b9654e08b8118bcbfead8e9dd2fdb9b053e9bdfb665285c78718f726d0b3d6c37e116428ec9ac9db2637259e4e8d6402bbada46c6bdb03985e19a82e9b4e57de1b025a3cb1f850beae7e8da9941655825bce0e89d536b6ee9064865b1a85c185e9fc9cb7f435de13d44773c00eed442a286e4ab807e3cab4dc3441d1b7d2af693812ae8b39652bb8c835fc895d13d6da93541afeadeee450475c29f3b2dfa8ef1c1e2547463b2cc2f0ff7a42ac4dd35e25c4fa030d2d2766fbe9f2d04c1304671747bace2f7dd55142bfa60f8cbc968bfc3d7a342152dc684a0fb5a32c0962a62b5220ac0f72add9d8b84d6cc76b97d03245e01fc8da3414a49bb4075d3488f29b56dc42ba69e3b58529448c943ecfd98b3784a39d0b8609a8fb945e757f4569f53bd2cf80f7f638acf5b67fe9c560a3b7b0cf7e0398f31aa8b03cf9c62b24296b6d8596b694469a02686c38daa16a1ef86e012d61a2f7de1693a5c00b3685175caec3c67146477eba54830f1d546cb18a553779aa46adb4f2010e33f3def847c7d89b51a8462b227605f6c920fd558a6daf64bc98682e508ae960c0c571870e603ba1fce0c13d53176f353fd319959e13db93eae1359f06e3dd4767c04f824cf34ec7bf8f60161ba1a615db82852eca9e3869afa711ab9a090660b0dc6cfbea310dda77e02310fbaeacd2636f975838c2dbcdbe9ac2cd85cee28f5e3f0c73abf62f9fa02cd79a7606b7ba855db68a07848b057c3aaf38f1a70086e14616f6f88305a1f9ce6b41378a620d4db3e0e7e1d421590dccaeff86212e232eeb5eb8a8d33a8c9b25ae88f3a7bd5032b4efa68f8af3186a02ffcbf5456f12beccace94c81c360cc4a0dcc642b59f991eec68c59af78139ca60b96d6a18e9535f8995e89bd2cf6a0aef3acffd33d1c0c1b79b66414a91d9f65b2b4ec65844b96f725d2b4b0c309f3eb9d714e9dd939bbdfd85ce8fb43679aeab13f6c29549949503c9466dbd337c4cdde46d6eacd15f21f4d8fdeaa627a47884c88a9c85f0b731d271a8ea7cb9e04a4a149c23c10f56b3a0476dc77a999d6e4f813e4b0f805e2a693e2ae4ae0ecc423c9ba5d17b42e691abf83784a582f2b1fd85d1e0a27ba38a500963568b2450363d2c5e3f7b8ba3e5b56e4e9f745a3a710bf2ae233c303068c532ce78ff031e6ab28b705dd94d7db4500909edb5626b8c9bd5ff4f0b4741388f0b91563ee516934c013e901572cba005ac5c535f4f107903be9af7b2793dfb61b5070facbe71eefe1b5600f975c8c38c3a2350d78beadfecb78e981164ae8bc866e732972d3ceef4aac68e15861f9b881d9b51b4edece150bc124b07645defb4202ef5d0e0962db98cae6ed459561c93c74c20bd64362e4f4fffc389a6cd80514604ff22eecc10c9cbc7981d19a8102b24146354c463107c9dc070e29e70df3578022acf72289ef071ab9f9402a544d0399f1b1e5f206b6d46d445f6d612a490e72918e00c853eda8493bef511149e80c9ab56e8b4b8cba3987249f77d060e61760e5792ac321c987c03c2606e9393a7970212992cdbd16448078d5039d4c2c3199714f53278f4f7b1d2e514cf95bdfc078b8bb0db659cb2c3f5cc02890ea84f05d414c88d2db9e9f8455659b9fa6254405317245fa070d6970cafb4dadb2522b490a5c8e02fe973a8cdbfbfbdbfb01535099ffba3d3896bc4d1189fc570c3e6fdc6469265b8da912772e75dd62ab71be507f700d56cac5e68fd6b57ec166168ab5258a69625c142a5b1b3519f94be1bde5e51d3bd8ea0c12d5af2fe4615b1b7bd4a96628a4fabc65925ff09718f63bbebaad98f89bd9543a27b3ff3b5d8bfa89f941a5eb8cc005ccd4a705190e1c9dc6a9f4264e5ee658520a4438e92de854bffc39f8dc7dfbb5de4f14ba63ea16a37d14a7b4610f95b6cffd55e4679b29cedbdf20e7bd16da822fad910c359ee3a68e48aae6e769b0e291d5d3aa3e2ca9d8d23abe8a1d5349f4991e9300852cc0befb20c2fc0d169306b260763344024f8092cbcc24c6807363e9fc548a30d5faab3a94b2af0782a2942be80c45d8b0587efd587394ef33c33022436e285806ddffdd32fe36345c3c38ed8d680abeb7a028b44ee6f94d060a14c7019bb6af1f1b5f0a562957d19826d8cc216f9b908c989ccd5415e3525dfe9422ffb5b50b7cc3083dc325544751e5683535d7439d3da2b0bb73bea551dd99e04e0e793804f4774eb6b1daf781d9caa5128274e599e847862fe309027813d3e4eda0bbeb7201856a5c5d8370e44dabff0bb229c723ba0a6bcf29c44536147de11b7835991018100105bd4329217f7386903fe8e7363cd7b3e893244e245e0a187467664c05b0be1fd429722b9b9a5e3198147fad72776e8a63aab9054fa9d259af0198d088d71d132e6068676a8e9ebb0f616b51ee34aac39c2c2221c71124017270d75ff4a048363c389e04e9b440ad2032a381ac2cfc54f409caa791e65ee4f5d6cd035008f219b88a803a7382ae447bf65a3df2176b25b3b7b67dabe34decd9a1384dc7a003916ca8fbcb29b3ad6fd8eac5bbbaa3bdfa6c6a3ad9427c4f3ed79fea26e14c8ce5fa3b4f82c5f7b6d2125916753a7b92ce9b46d45 diff --git a/numpy/random/src/dsfmt/96-bit-jump.poly.txt b/numpy/random/src/dsfmt/96-bit-jump.poly.txt deleted file mode 100644 index 15c68d155..000000000 --- a/numpy/random/src/dsfmt/96-bit-jump.poly.txt +++ /dev/null @@ -1,2 +0,0 @@ -jump polynomial: -288da521f7244e5f62bf26692bdd1fcdfd38a850addb02d98bd358367cb78c71348f3e163522e0e30e4feb90aa210dd793c94d151fa89aa319911aa42b511ea568a503d595c6d9bcd37317a37ef5679700a5b67f29df72451770fc1eb8c97427cdd9825c23f32dcd5c4fb117a4f5982a3bee8f16595d021165cd9688db342e360e222c9855c8306fd7b5fc82e62e3b1765e7f3319da9da66c325b030bd6175876efc70636306cd2de31a299ca20e9eb1d5063bcbff0ba282aff4737a5b1585cd940ae9cd45fda222308341e0d5588e81b42c4e0574deeb2d80b84c00cb3f8a7ae6278462e1994b83a25b33aa0dc74d5d3d057dabfd6a8a82d7dfb6bb66a223bc46dca2b5fb1885e6ab80fddcd6578b32c21c4a9d761cb9798800c921d56ee356c491454e15956e68ef566c1706fcdfb6a6828ec1fb93db455525828e8741a371c14959658b99bbd358a162230ee451a8432df279e4ba0d3a493954359a5390b16475540578270b168054fefb1e20948d4d20c88005ed09e671b6a94b8ea929f72e7b2f85af4098a92d742b64964ea6b7773b2c20b22a0ff35bd9367c3160b549411230e15a987f361e04daac49d2fe4c7c9371d84bf270d8f33a62680b2ee014bf5be691aa0d82e66e885eaa832a241aff8a09c435ac2b0698bc3865c5125d498a4ffadd31d5f2c6aee6496fdc6c13055b53e7f65a683ef539b6e8ea6e21429a11ff74ccef09ee83eac1b5ddaf1b77fed786fd20e8cbb3e8877b7f80a24fef31a9f8d8108099c50abc886f7ab20c4396bf51af1b806003185eaf2061f55036e26f498b92aabadfb6b5bed2d116b32ae387a692058e6160a9508dc44c971454d9a357ba115278861be0aeaa0926d773c5d4563e37dffcfed8bbf44b79e384650b7eff30aae73154a2ef130cee1eaf32d944e5472ae217789c03276deb8290c73dd5cde0b6dce9b28cbb73776e3e52260768f37a38db3d1c26c3c065b641c7a825edc155d947118d8b6ff8c8818440088724261ca83fa866aa4159fbffac8c28c8a7ca5f1e2fde74b4908c8215cbde20190bdf0de1d5a05a2c116a66eeadcafd643098e74ec3e82b3658c0c4fd7c7797d567b8db3d6b67ca262d713dbf45cc80b89f279be0991f276a4501d2ea6222026caa7e6fbcf4d02fdf65d8f439f88cfb1121d1b0f4dd670d146b93e099a32869448d837e347229745e5c30f1521b0c477b2062c9c8f631dcd29664eec7f28bdcac2a1ca2deabbbc89b21824ba92a61eeb4c5dd62b20c346134d842bcfc189f0e2244bfb8596c38c88d5bd4447fcd30890a4acf71801a6bcf5d806073b9ca730db53e160a79516333748708dd5e5be37d75e90e64d14ddf5ccc9390ae67cbba39916ce9b3b6b1d19378e4bd36ef0c9e62570394051cc273122e0790e381b20e083beca6e88bc2fa8bde22c89689b4b710c775cd9065158b48bf855fc3a3e80693908046ea1da9c558f024ea5ea69d65f40317fc3c8bab7606bf7edf17fcaeb706756c4de38319a51fc24c80a2baccef40f4354f5147fb91c9b1b91011d146da7eeb426d25de17d7e39ee53735ef28b64d5e6e5444676f679a8e4e91314377c7606531b90dc1cd8cf2d929ed9137fdc0d6b6db62e764ef200a439d742b51425b818231ed799a53768119466cc42296ce36f52507411709cd9d622c1f899291db27104589a5e4509d9818cef983d6a13ce1247184c4551c38bd319aa068cd9c0b96f9e08b4a7bd7852c129d4126676cbcb30ae029b0e2cec3c58c110fecae7967fca0752874e2fcc077084dd4d87b4dee632b621f35cb9f22639ab4de86089c56cabb1aa8a4fedbc5d8673cca12b4479dca5dc6f3048bb654fd59c464b033960177a4d6c2ee38b124d193cd5da3cbc5aa5ebdf7197f5d3f3008f01effb599b6d238044f2ee204bf372685256813278ca7a7d145445f5d5eb5490d99ee89d587fe446e1600d7daf9553b7d709bd796c59757e2f7c26cb71990a51ffc49120f358ef4927729b7e8c2cf9ad35a8017667625a39e2d12c8b0dd5de3d3797d9167ac9302a122e0f066a2552a7278a632b043089a1f5a127ce0bc8e688c545f20bca467fd58d9f9b8234f2be7a28152ab4d39ba8d6c4471d574571aa1a4a3aca16f46ac74e664899cf24502334aec637a9c392ba1618340080bfaed493abaa2f684ffb2bc1da86d804305772e30893f67235d5ce8180ef60490e46f131465d249a39df1aaed30621169c0f7db0b6a6adcab975ec00ca68b2fc06f443cfb8933b326f106670e0e528f22ff96d06b2a6d917258d29a87cf637f37c3661b42462c11c73c1cd0a6774e4e10157b02c0cfd9148ad8a23775653db3dec41fd6d1d5eb54a3e5e144a489a18412f21eb96a72b6751482882d636a1cd0057ea48d1a985472b0c301c4a5b24b0152bdc50a47126a01a1d0665270cdbf2ed852e0f61515bad6e58973329d269bf057ffa52331dde4700b926f635cdf49234f4aaabbabca5386526568fc3a80b1d8a9e1566e21bf889546379289263d978e95de390c4bbdb5a813535b7186f661f5c283400adb3bcf9365d45abb6088b729a98789265f1e8e684e2b2c15b4ce53230abc5e5bf6895827c95842e0849fc7fe56b7c65f075baded0f2e172c8e1088219615b8697ff4a3c6f5f708e498e6c7312680f214a877061511d3b85087545fc26708d039e977a0157b95ceba40497cc2dd1b1b1394eb98f651d701dfbc3583c159da8bae76a25db213211d6191d83c5d7a3d4b2320b8a305b5d90e04d6e3de161fff5dae7e3d4f6c9e1cf55cb5ac1677d9cfd560db2209be2b9566a267aad9cf4b9e03c221ba5b0f02cabb50cef19180ba329691d114da670d7b85e36c0b6b7d81613bd31350dfe225050861e90043ac2334478f52584a1a8809f76d40af36da02549c54da164487f0b7f823cff797db1c90f2f1c431eca97fc8bcdfb44c30cd1643c893d5e33aa56cbc0a0c38f5c8f6bb37483d13b85b659ac51ce05311bba19c97772e81c2315a009e80f591c82445d493dc3b5fb12c52e8c50c6260200b0d77092bf19aabce491b2c511fca2a4ebd99b446e55f453314297723894f0b223868ef80f5964468afa3a5e6aa41f2128e6de893bba2cbde9bea91ba97bda18a01905ce8d2e85e6011cc0550f5ae9121361fdec1ec97a6e6892a68732a69147476d54babaa564b883baad7100eb1092a1aa28a29f67e6b53deab511e53eada87cfe1bb6e3c6a215fd8d729840a5b5ac452cfd8acb9be7818d2c806c3e0cedd0847ddf9a5906bf1a0fa4da44ccea829a1a5350d5db0241a221b97e5cd5ba15e58b402529317c854fbcda86a562d4ee44a34193513647af0a3bc9f90ababc1dbbfd9aba8d3dcc39463473ca6bc0e1dc736ba712eef42dee80e31e7d8abe23f98e91ab875d0553bc24be9cb1d9484812c0b038cb177ad52064328e17f8ca7c8737902d964017e3aaae66161270dac21de42a6f60d10d89c1556916a249a130752bb7c7783b93a59d9f5456745ecc512f497b5a31be2678b9587628cb45dae2f5f6bde7ea4500c1ba961e2 diff --git a/numpy/random/src/dsfmt/LICENSE.md b/numpy/random/src/dsfmt/LICENSE.md deleted file mode 100644 index d59568f6b..000000000 --- a/numpy/random/src/dsfmt/LICENSE.md +++ /dev/null @@ -1,34 +0,0 @@ -# DSFMT - -Copyright (c) 2007, 2008, 2009 Mutsuo Saito, Makoto Matsumoto -and Hiroshima University. -Copyright (c) 2011, 2002 Mutsuo Saito, Makoto Matsumoto, Hiroshima -University and The University of Tokyo. -All rights reserved. - -Redistribution and use in source and binary forms, with or without -modification, are permitted provided that the following conditions are -met: - -* Redistributions of source code must retain the above copyright - notice, this list of conditions and the following disclaimer. -* 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. -* Neither the name of the Hiroshima University nor the names of - its contributors may 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. diff --git a/numpy/random/src/dsfmt/calc-jump.cpp b/numpy/random/src/dsfmt/calc-jump.cpp deleted file mode 100644 index 495b2797c..000000000 --- a/numpy/random/src/dsfmt/calc-jump.cpp +++ /dev/null @@ -1,81 +0,0 @@ -/** - * @file calc-jump.cpp - * - * @brief calc jump function. - * - * @author Mutsuo Saito (Hiroshima University) - * @author Makoto Matsumoto (The University of Tokyo) - * - * Copyright (C) 2012 Mutsuo Saito, Makoto Matsumoto, - * Hiroshima University and The University of Tokyo. - * All rights reserved. - * - * The 3-clause BSD License is applied to this software, see - * LICENSE.txt - * - * Compile: - * g++ calc-jump.cpp -o calc-jump -lntl - * - * Compute polynomial for 2^128 steps: - * ./calc-jump 340282366920938463463374607431768211456 poly.19937.txt - * - */ -#include <iostream> -#include <fstream> -#include <iomanip> -#include <sstream> -#include <string> -#include <inttypes.h> -#include <stdint.h> -#include <time.h> -#include <NTL/GF2X.h> -#include <NTL/vec_GF2.h> -#include <NTL/ZZ.h> -#include "dSFMT-calc-jump.hpp" - -using namespace NTL; -using namespace std; -using namespace dsfmt; - -static void read_file(GF2X& lcmpoly, long line_no, const string& file); - -int main(int argc, char * argv[]) { - if (argc <= 2) { - cout << argv[0] << " jump-step poly-file" << endl; - cout << " jump-step: a number between zero and 2^{DSFMT_MEXP}-1.\n" - << " large decimal number is allowed." << endl; - cout << " poly-file: one of poly.{MEXP}.txt " - << "file" << endl; - return -1; - } - string step_string = argv[1]; - string filename = argv[2]; - long no = 0; - GF2X lcmpoly; - read_file(lcmpoly, no, filename); - ZZ step; - stringstream ss(step_string); - ss >> step; - string jump_str; - calc_jump(jump_str, step, lcmpoly); - cout << "jump polynomial:" << endl; - cout << jump_str << endl; - return 0; -} - - -static void read_file(GF2X& lcmpoly, long line_no, const string& file) -{ - ifstream ifs(file.c_str()); - string line; - for (int i = 0; i < line_no; i++) { - ifs >> line; - ifs >> line; - } - if (ifs) { - ifs >> line; - line = ""; - ifs >> line; - } - stringtopoly(lcmpoly, line); -} diff --git a/numpy/random/src/dsfmt/dSFMT-benchmark.c b/numpy/random/src/dsfmt/dSFMT-benchmark.c deleted file mode 100644 index af29d0e1f..000000000 --- a/numpy/random/src/dsfmt/dSFMT-benchmark.c +++ /dev/null @@ -1,43 +0,0 @@ -/* - * - * cl dsfmt-benchmark.c dSFMT.c /Ox -DHAVE_SSE2 - * - * gcc dSFMT-benchmark.c dSFMT.c -O3 -DHAVE_SSE2 -DDSFMT_MEXP=19937 -o - * dSFMT-benchmark - */ -#include <inttypes.h> -#include <time.h> - -#include "dSFMT.h" - - -#define N 1000000000 - -int main() { - int i, j; - uint32_t seed = 0xDEADBEAF; - uint64_t count = 0, sum = 0; - dsfmt_t state; - double buffer[DSFMT_N64]; - - uint64_t out; - uint64_t *tmp; - dsfmt_init_gen_rand(&state, seed); - clock_t begin = clock(); - for (i = 0; i < N / (DSFMT_N64 / 2); i++) { - dsfmt_fill_array_close_open(&state, &buffer[0], DSFMT_N64); - for (j = 0; j < DSFMT_N64; j += 2) { - tmp = (uint64_t *)&buffer[j]; - out = (*tmp >> 16) << 32; - tmp = (uint64_t *)&buffer[j + 1]; - out |= (*tmp >> 16) & 0xffffffff; - sum += out; - count++; - } - } - clock_t end = clock(); - double time_spent = (double)(end - begin) / CLOCKS_PER_SEC; - printf("0x%" PRIx64 "\ncount: %" PRIu64 "\n", sum, count); - printf("%" PRIu64 " randoms per second\n", - (uint64_t)(N / time_spent) / 1000000 * 1000000); -} diff --git a/numpy/random/src/dsfmt/dSFMT-calc-jump.hpp b/numpy/random/src/dsfmt/dSFMT-calc-jump.hpp deleted file mode 100644 index b960826be..000000000 --- a/numpy/random/src/dsfmt/dSFMT-calc-jump.hpp +++ /dev/null @@ -1,106 +0,0 @@ -#pragma once -#ifndef DSFMT_CALC_JUMP_HPP -#define DSFMT_CALC_JUMP_HPP -/** - * @file dSFMT-calc-jump.hpp - * - * @brief functions for calculating jump polynomial. - * - * @author Mutsuo Saito (Hiroshima University) - * @author Makoto Matsumoto (The University of Tokyo) - * - * Copyright (C) 2012 Mutsuo Saito, Makoto Matsumoto, - * Hiroshima University and The University of Tokyo. - * All rights reserved. - * - * The 3-clause BSD License is applied to this software, see - * LICENSE.txt - */ -#include <iostream> -#include <iomanip> -#include <sstream> -#include <NTL/GF2X.h> - -namespace dsfmt { -/** - * converts polynomial to string for convenient use in C language. - * @param x output string - * @param polynomial input polynomial - */ - static inline void polytostring(std::string& x, NTL::GF2X& polynomial) - { - using namespace NTL; - using namespace std; - - long degree = deg(polynomial); - int buff; - stringstream ss; - for (int i = 0; i <= degree; i+=4) { - buff = 0; - for (int j = 0; j < 4; j++) { - if (IsOne(coeff(polynomial, i + j))) { - buff |= 1 << j; - } else { - buff &= (0x0f ^ (1 << j)); - } - } - ss << hex << buff; - } - ss << flush; - x = ss.str(); - } - -/** - * converts string to polynomial - * @param str string - * @param poly output polynomial - */ - static inline void stringtopoly(NTL::GF2X& poly, std::string& str) - { - using namespace NTL; - using namespace std; - - stringstream ss(str); - char c; - long p = 0; - clear(poly); - while(ss) { - ss >> c; - if (!ss) { - break; - } - if (c >= 'a') { - c = c - 'a' + 10; - } else { - c = c - '0'; - } - for (int j = 0; j < 4; j++) { - if (c & (1 << j)) { - SetCoeff(poly, p, 1); - } else { - SetCoeff(poly, p, 0); - } - p++; - } - } - } - -/** - * calculate the jump polynomial. - * SFMT generates 4 32-bit integers from one internal state. - * @param jump_str output string which represents jump polynomial. - * @param step jump step of internal state - * @param characteristic polynomial - */ - static inline void calc_jump(std::string& jump_str, - NTL::ZZ& step, - NTL::GF2X& characteristic) - { - using namespace NTL; - using namespace std; - GF2X jump; - PowerXMod(jump, step, characteristic); - polytostring(jump_str, jump); - } -} -#endif diff --git a/numpy/random/src/dsfmt/dSFMT-common.h b/numpy/random/src/dsfmt/dSFMT-common.h deleted file mode 100644 index 30c26c08b..000000000 --- a/numpy/random/src/dsfmt/dSFMT-common.h +++ /dev/null @@ -1,115 +0,0 @@ -#pragma once -/** - * @file dSFMT-common.h - * - * @brief SIMD oriented Fast Mersenne Twister(SFMT) pseudorandom - * number generator with jump function. This file includes common functions - * used in random number generation and jump. - * - * @author Mutsuo Saito (Hiroshima University) - * @author Makoto Matsumoto (The University of Tokyo) - * - * Copyright (C) 2006, 2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima - * University. - * Copyright (C) 2012 Mutsuo Saito, Makoto Matsumoto, Hiroshima - * University and The University of Tokyo. - * All rights reserved. - * - * The 3-clause BSD License is applied to this software, see - * LICENSE.txt - */ -#ifndef DSFMT_COMMON_H -#define DSFMT_COMMON_H - -#include "dSFMT.h" - -#if defined(HAVE_SSE2) -# include <emmintrin.h> -union X128I_T { - uint64_t u[2]; - __m128i i128; -}; -union X128D_T { - double d[2]; - __m128d d128; -}; -/** mask data for sse2 */ -static const union X128I_T sse2_param_mask = {{DSFMT_MSK1, DSFMT_MSK2}}; -#endif - -#if defined(HAVE_ALTIVEC) -inline static void do_recursion(w128_t *r, w128_t *a, w128_t * b, - w128_t *lung) { - const vector unsigned char sl1 = ALTI_SL1; - const vector unsigned char sl1_perm = ALTI_SL1_PERM; - const vector unsigned int sl1_msk = ALTI_SL1_MSK; - const vector unsigned char sr1 = ALTI_SR; - const vector unsigned char sr1_perm = ALTI_SR_PERM; - const vector unsigned int sr1_msk = ALTI_SR_MSK; - const vector unsigned char perm = ALTI_PERM; - const vector unsigned int msk1 = ALTI_MSK; - vector unsigned int w, x, y, z; - - z = a->s; - w = lung->s; - x = vec_perm(w, (vector unsigned int)perm, perm); - y = vec_perm(z, (vector unsigned int)sl1_perm, sl1_perm); - y = vec_sll(y, sl1); - y = vec_and(y, sl1_msk); - w = vec_xor(x, b->s); - w = vec_xor(w, y); - x = vec_perm(w, (vector unsigned int)sr1_perm, sr1_perm); - x = vec_srl(x, sr1); - x = vec_and(x, sr1_msk); - y = vec_and(w, msk1); - z = vec_xor(z, y); - r->s = vec_xor(z, x); - lung->s = w; -} -#elif defined(HAVE_SSE2) -/** - * This function represents the recursion formula. - * @param r output 128-bit - * @param a a 128-bit part of the internal state array - * @param b a 128-bit part of the internal state array - * @param d a 128-bit part of the internal state array (I/O) - */ -inline static void do_recursion(w128_t *r, w128_t *a, w128_t *b, w128_t *u) { - __m128i v, w, x, y, z; - - x = a->si; - z = _mm_slli_epi64(x, DSFMT_SL1); - y = _mm_shuffle_epi32(u->si, SSE2_SHUFF); - z = _mm_xor_si128(z, b->si); - y = _mm_xor_si128(y, z); - - v = _mm_srli_epi64(y, DSFMT_SR); - w = _mm_and_si128(y, sse2_param_mask.i128); - v = _mm_xor_si128(v, x); - v = _mm_xor_si128(v, w); - r->si = v; - u->si = y; -} -#else -/** - * This function represents the recursion formula. - * @param r output 128-bit - * @param a a 128-bit part of the internal state array - * @param b a 128-bit part of the internal state array - * @param lung a 128-bit part of the internal state array (I/O) - */ -inline static void do_recursion(w128_t *r, w128_t *a, w128_t * b, - w128_t *lung) { - uint64_t t0, t1, L0, L1; - - t0 = a->u[0]; - t1 = a->u[1]; - L0 = lung->u[0]; - L1 = lung->u[1]; - lung->u[0] = (t0 << DSFMT_SL1) ^ (L1 >> 32) ^ (L1 << 32) ^ b->u[0]; - lung->u[1] = (t1 << DSFMT_SL1) ^ (L0 >> 32) ^ (L0 << 32) ^ b->u[1]; - r->u[0] = (lung->u[0] >> DSFMT_SR) ^ (lung->u[0] & DSFMT_MSK1) ^ t0; - r->u[1] = (lung->u[1] >> DSFMT_SR) ^ (lung->u[1] & DSFMT_MSK2) ^ t1; -} -#endif -#endif diff --git a/numpy/random/src/dsfmt/dSFMT-jump.c b/numpy/random/src/dsfmt/dSFMT-jump.c deleted file mode 100644 index 1832bb885..000000000 --- a/numpy/random/src/dsfmt/dSFMT-jump.c +++ /dev/null @@ -1,184 +0,0 @@ -/** - * @file dSFMT-jump.c - * - * @brief do jump using jump polynomial. - * - * @author Mutsuo Saito (Hiroshima University) - * @author Makoto Matsumoto (The University of Tokyo) - * - * Copyright (C) 2012 Mutsuo Saito, Makoto Matsumoto, - * Hiroshima University and The University of Tokyo. - * All rights reserved. - * - * The 3-clause BSD License is applied to this software, see - * LICENSE.txt - */ - -#include <assert.h> -#include <stdlib.h> -#include <ctype.h> -#include <string.h> -#include "dSFMT-params.h" -#include "dSFMT.h" -#include "dSFMT-jump.h" -#include "dSFMT-common.h" - -#if defined(__cplusplus) -extern "C" { -#endif - - struct FIX_T { - int mexp; - uint64_t fix[4]; - }; - - struct FIX_T fix_table[] = { - {521, {UINT64_C(0x3fff56977f035125), - UINT64_C(0x3ff553857b015035), - UINT64_C(0x4034434434434434), - UINT64_C(0x0140151151351371)}}, - {1279, {UINT64_C(0x3ff87befce70e89f), - UINT64_C(0x3ff5f6afa3c60868), - UINT64_C(0xa4ca4caccaccacdb), - UINT64_C(0x40444444444c44c4)}}, - {4253, {UINT64_C(0x3ff85a66da51a81a), - UINT64_C(0x3ff4f4aeab9688eb), - UINT64_C(0x20524524534d34d3), - UINT64_C(0xc9cc9cc9cc9ccdcf)}}, - {216091, {UINT64_C(0x3ff096d54a871071), - UINT64_C(0x3ffafa9bfbd5d55d), - UINT64_C(0x0470470470573573), - UINT64_C(0x0250250251259259)}}, - {0} - }; - - inline static void next_state(dsfmt_t * dsfmt); - -#if defined(HAVE_SSE2) -/** - * add internal state of src to dest as F2-vector. - * @param dest destination state - * @param src source state - */ - inline static void add(dsfmt_t *dest, dsfmt_t *src) { - int dp = dest->idx / 2; - int sp = src->idx / 2; - int diff = (sp - dp + DSFMT_N) % DSFMT_N; - int p; - int i; - for (i = 0; i < DSFMT_N - diff; i++) { - p = i + diff; - dest->status[i].si - = _mm_xor_si128(dest->status[i].si, src->status[p].si); - } - for (i = DSFMT_N - diff; i < DSFMT_N; i++) { - p = i + diff - DSFMT_N; - dest->status[i].si - = _mm_xor_si128(dest->status[i].si, src->status[p].si); - } - dest->status[DSFMT_N].si - = _mm_xor_si128(dest->status[DSFMT_N].si, - src->status[DSFMT_N].si); - } -#else - inline static void add(dsfmt_t *dest, dsfmt_t *src) { - int dp = dest->idx / 2; - int sp = src->idx / 2; - int diff = (sp - dp + DSFMT_N) % DSFMT_N; - int p; - int i; - for (i = 0; i < DSFMT_N - diff; i++) { - p = i + diff; - dest->status[i].u[0] ^= src->status[p].u[0]; - dest->status[i].u[1] ^= src->status[p].u[1]; - } - for (; i < DSFMT_N; i++) { - p = i + diff - DSFMT_N; - dest->status[i].u[0] ^= src->status[p].u[0]; - dest->status[i].u[1] ^= src->status[p].u[1]; - } - dest->status[DSFMT_N].u[0] ^= src->status[DSFMT_N].u[0]; - dest->status[DSFMT_N].u[1] ^= src->status[DSFMT_N].u[1]; - } -#endif - -/** - * calculate next state - * @param dsfmt dSFMT internal state - */ - inline static void next_state(dsfmt_t * dsfmt) { - int idx = (dsfmt->idx / 2) % DSFMT_N; - w128_t * lung; - w128_t * pstate = &dsfmt->status[0]; - - lung = &pstate[DSFMT_N]; - do_recursion(&pstate[idx], - &pstate[idx], - &pstate[(idx + DSFMT_POS1) % DSFMT_N], - lung); - dsfmt->idx = (dsfmt->idx + 2) % DSFMT_N64; - } - - inline static void add_fix(dsfmt_t * dsfmt) { - int i; - int index = -1; - for (i = 0; fix_table[i].mexp != 0; i++) { - if (fix_table[i].mexp == DSFMT_MEXP) { - index = i; - } - if (fix_table[i].mexp > DSFMT_MEXP) { - break; - } - } - if (index < 0) { - return; - } - for (i = 0; i < DSFMT_N; i++) { - dsfmt->status[i].u[0] ^= fix_table[index].fix[0]; - dsfmt->status[i].u[1] ^= fix_table[index].fix[1]; - } - dsfmt->status[DSFMT_N].u[0] ^= fix_table[index].fix[2]; - dsfmt->status[DSFMT_N].u[1] ^= fix_table[index].fix[3]; - } - -/** - * jump ahead using jump_string - * @param dsfmt dSFMT internal state input and output. - * @param jump_string string which represents jump polynomial. - */ - void dSFMT_jump(dsfmt_t * dsfmt, const char * jump_string) { - dsfmt_t work; - int index = dsfmt->idx; - int bits; - int i; - int j; - memset(&work, 0, sizeof(dsfmt_t)); - add_fix(dsfmt); - dsfmt->idx = DSFMT_N64; - - for (i = 0; jump_string[i] != '\0'; i++) { - bits = jump_string[i]; - assert(isxdigit(bits)); - bits = tolower(bits); - if (bits >= 'a' && bits <= 'f') { - bits = bits - 'a' + 10; - } else { - bits = bits - '0'; - } - bits = bits & 0x0f; - for (j = 0; j < 4; j++) { - if ((bits & 1) != 0) { - add(&work, dsfmt); - } - next_state(dsfmt); - bits = bits >> 1; - } - } - *dsfmt = work; - add_fix(dsfmt); - dsfmt->idx = index; - } - -#if defined(__cplusplus) -} -#endif diff --git a/numpy/random/src/dsfmt/dSFMT-jump.h b/numpy/random/src/dsfmt/dSFMT-jump.h deleted file mode 100644 index 689f9499a..000000000 --- a/numpy/random/src/dsfmt/dSFMT-jump.h +++ /dev/null @@ -1,29 +0,0 @@ -#pragma once -#ifndef DSFMT_JUMP_H -#define DSFMT_JUMP_H -/** - * @file SFMT-jump.h - * - * @brief jump header file. - * - * @author Mutsuo Saito (Hiroshima University) - * @author Makoto Matsumoto (The University of Tokyo) - * - * Copyright (C) 2012 Mutsuo Saito, Makoto Matsumoto, - * Hiroshima University and The University of Tokyo. - * All rights reserved. - * - * The 3-clause BSD License is applied to this software, see - * LICENSE.txt - */ -#if defined(__cplusplus) -extern "C" { -#endif - -#include "dSFMT.h" -void dSFMT_jump(dsfmt_t *dsfmt, const char *jump_str); - -#if defined(__cplusplus) -} -#endif -#endif diff --git a/numpy/random/src/dsfmt/dSFMT-params.h b/numpy/random/src/dsfmt/dSFMT-params.h deleted file mode 100644 index aa0247800..000000000 --- a/numpy/random/src/dsfmt/dSFMT-params.h +++ /dev/null @@ -1,87 +0,0 @@ -#ifndef DSFMT_PARAMS_H -#define DSFMT_PARAMS_H - -#include "dSFMT.h" - -/*---------------------- - the parameters of DSFMT - following definitions are in dSFMT-paramsXXXX.h file. - ----------------------*/ -/** the pick up position of the array. -#define DSFMT_POS1 122 -*/ - -/** the parameter of shift left as four 32-bit registers. -#define DSFMT_SL1 18 - */ - -/** the parameter of shift right as four 32-bit registers. -#define DSFMT_SR1 12 -*/ - -/** A bitmask, used in the recursion. These parameters are introduced - * to break symmetry of SIMD. -#define DSFMT_MSK1 (uint64_t)0xdfffffefULL -#define DSFMT_MSK2 (uint64_t)0xddfecb7fULL -*/ - -/** These definitions are part of a 128-bit period certification vector. -#define DSFMT_PCV1 UINT64_C(0x00000001) -#define DSFMT_PCV2 UINT64_C(0x00000000) -*/ - -#define DSFMT_LOW_MASK UINT64_C(0x000FFFFFFFFFFFFF) -#define DSFMT_HIGH_CONST UINT64_C(0x3FF0000000000000) -#define DSFMT_SR 12 - -/* for sse2 */ -#if defined(HAVE_SSE2) - #define SSE2_SHUFF 0x1b -#elif defined(HAVE_ALTIVEC) - #if defined(__APPLE__) /* For OSX */ - #define ALTI_SR (vector unsigned char)(4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4) - #define ALTI_SR_PERM \ - (vector unsigned char)(15,0,1,2,3,4,5,6,15,8,9,10,11,12,13,14) - #define ALTI_SR_MSK \ - (vector unsigned int)(0x000fffffU,0xffffffffU,0x000fffffU,0xffffffffU) - #define ALTI_PERM \ - (vector unsigned char)(12,13,14,15,8,9,10,11,4,5,6,7,0,1,2,3) - #else - #define ALTI_SR {4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4} - #define ALTI_SR_PERM {15,0,1,2,3,4,5,6,15,8,9,10,11,12,13,14} - #define ALTI_SR_MSK {0x000fffffU,0xffffffffU,0x000fffffU,0xffffffffU} - #define ALTI_PERM {12,13,14,15,8,9,10,11,4,5,6,7,0,1,2,3} - #endif -#endif - -#if DSFMT_MEXP == 521 - #include "dSFMT-params521.h" -#elif DSFMT_MEXP == 1279 - #include "dSFMT-params1279.h" -#elif DSFMT_MEXP == 2203 - #include "dSFMT-params2203.h" -#elif DSFMT_MEXP == 4253 - #include "dSFMT-params4253.h" -#elif DSFMT_MEXP == 11213 - #include "dSFMT-params11213.h" -#elif DSFMT_MEXP == 19937 - #include "dSFMT-params19937.h" -#elif DSFMT_MEXP == 44497 - #include "dSFMT-params44497.h" -#elif DSFMT_MEXP == 86243 - #include "dSFMT-params86243.h" -#elif DSFMT_MEXP == 132049 - #include "dSFMT-params132049.h" -#elif DSFMT_MEXP == 216091 - #include "dSFMT-params216091.h" -#else -#ifdef __GNUC__ - #error "DSFMT_MEXP is not valid." - #undef DSFMT_MEXP -#else - #undef DSFMT_MEXP -#endif - -#endif - -#endif /* DSFMT_PARAMS_H */ diff --git a/numpy/random/src/dsfmt/dSFMT-params19937.h b/numpy/random/src/dsfmt/dSFMT-params19937.h deleted file mode 100644 index a600b0dbc..000000000 --- a/numpy/random/src/dsfmt/dSFMT-params19937.h +++ /dev/null @@ -1,40 +0,0 @@ -#ifndef DSFMT_PARAMS19937_H -#define DSFMT_PARAMS19937_H - -/* #define DSFMT_N 191 */ -/* #define DSFMT_MAXDEGREE 19992 */ -#define DSFMT_POS1 117 -#define DSFMT_SL1 19 -#define DSFMT_MSK1 UINT64_C(0x000ffafffffffb3f) -#define DSFMT_MSK2 UINT64_C(0x000ffdfffc90fffd) -#define DSFMT_MSK32_1 0x000ffaffU -#define DSFMT_MSK32_2 0xfffffb3fU -#define DSFMT_MSK32_3 0x000ffdffU -#define DSFMT_MSK32_4 0xfc90fffdU -#define DSFMT_FIX1 UINT64_C(0x90014964b32f4329) -#define DSFMT_FIX2 UINT64_C(0x3b8d12ac548a7c7a) -#define DSFMT_PCV1 UINT64_C(0x3d84e1ac0dc82880) -#define DSFMT_PCV2 UINT64_C(0x0000000000000001) -#define DSFMT_IDSTR "dSFMT2-19937:117-19:ffafffffffb3f-ffdfffc90fffd" - - -/* PARAMETERS FOR ALTIVEC */ -#if defined(__APPLE__) /* For OSX */ - #define ALTI_SL1 (vector unsigned char)(3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3) - #define ALTI_SL1_PERM \ - (vector unsigned char)(2,3,4,5,6,7,30,30,10,11,12,13,14,15,0,1) - #define ALTI_SL1_MSK \ - (vector unsigned int)(0xffffffffU,0xfff80000U,0xffffffffU,0xfff80000U) - #define ALTI_MSK (vector unsigned int)(DSFMT_MSK32_1, \ - DSFMT_MSK32_2, DSFMT_MSK32_3, DSFMT_MSK32_4) -#else /* For OTHER OSs(Linux?) */ - #define ALTI_SL1 {3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3} - #define ALTI_SL1_PERM \ - {2,3,4,5,6,7,30,30,10,11,12,13,14,15,0,1} - #define ALTI_SL1_MSK \ - {0xffffffffU,0xfff80000U,0xffffffffU,0xfff80000U} - #define ALTI_MSK \ - {DSFMT_MSK32_1, DSFMT_MSK32_2, DSFMT_MSK32_3, DSFMT_MSK32_4} -#endif - -#endif /* DSFMT_PARAMS19937_H */ diff --git a/numpy/random/src/dsfmt/dSFMT-poly.h b/numpy/random/src/dsfmt/dSFMT-poly.h deleted file mode 100644 index f8e15c3eb..000000000 --- a/numpy/random/src/dsfmt/dSFMT-poly.h +++ /dev/null @@ -1,53 +0,0 @@ -static const char * poly_128 = -"f4dfa6c62049d0776e0bf6f1e953f3aa38abb113df86be024eab3773ad5f2b82ead936022e656dff7e562691c59dd5f7d2" -"566b78d9669002503c4ddb1888a49f32333f515e6c60c4ecd221078ec6f26f0a90f4875067ca1f399a99775037adf90556" -"6e2c7e6b42131420f8f04f112c92621c9b1502f2a8aefad6c667904af62f0d55e02d396902d3b89450103c5ce5fe0408d9" -"7cbb864861b49e4e42048ff3310b48faac55095a7f422eea4aade752f947f947c6be0a0c665bdea099246ab9eff658ea8c" -"a468bf49d0227748367878de06d7bd86ea6708fcac6e252f5f00f04309b2aac3036b64afb39d990427c6c9f03477cc7e93" -"5c43c0e61bc161db8eb15516eee8cb377ecbc1849207990fb6778721b29bfe0d89bfda1b3772fa5b0b1f7ec3daf3605203" -"2285898c6f6396f55010c31f8201b7e2e51d94f920bfe57684c5415cc342cb39a0045d9793d13cf8646096daeb8bb9bfc2" -"0a90de8f2426da8733267a9b9674f32154e8f84a9932223a2ca3c787d0b66df6675febbdfcba2f9cef09c621c57e11098b" -"3289c77397aaae8b104642ffe0c4b75598efbc53745984d68b4d6656cae299ae2be55217a9a02b009ca7be32f47fbe434b" -"ce4914a34d0c9b0085bede9b8a99319c34660d66f0124b5a7714c4bf3cbfec3ee43ed817087168bad80133bebaeeb68cf7" -"929a24d1bb3de831a8340d220906ab04159cf94b21d5ee813bd7c80f10f01b43052af530917513b169254c25d6fcfe6cb4" -"20d6ce92f54886ef6eaf9a5ba35e893ff593834d05ddf28899e42d729c7df3d21ef036020789739366f0c11ec52ff92a0b" -"fd8ba69508e27b20fabb8217bd36b90e5aa918159ac87913bc7b46c04e366c23c92807fbe9c6a407e6a4db0b4fc23c3b6c" -"706b5ca058fe8c190f849f18d16d6b48b5ed760eb202fd566291a799420b9654e08b8118bcbfead8e9dd2fdb9b053e9bdf" -"b665285c78718f726d0b3d6c37e116428ec9ac9db2637259e4e8d6402bbada46c6bdb03985e19a82e9b4e57de1b025a3cb" -"1f850beae7e8da9941655825bce0e89d536b6ee9064865b1a85c185e9fc9cb7f435de13d44773c00eed442a286e4ab807e" -"3cab4dc3441d1b7d2af693812ae8b39652bb8c835fc895d13d6da93541afeadeee450475c29f3b2dfa8ef1c1e2547463b2" -"cc2f0ff7a42ac4dd35e25c4fa030d2d2766fbe9f2d04c1304671747bace2f7dd55142bfa60f8cbc968bfc3d7a342152dc6" -"84a0fb5a32c0962a62b5220ac0f72add9d8b84d6cc76b97d03245e01fc8da3414a49bb4075d3488f29b56dc42ba69e3b58" -"529448c943ecfd98b3784a39d0b8609a8fb945e757f4569f53bd2cf80f7f638acf5b67fe9c560a3b7b0cf7e0398f31aa8b" -"03cf9c62b24296b6d8596b694469a02686c38daa16a1ef86e012d61a2f7de1693a5c00b3685175caec3c67146477eba548" -"30f1d546cb18a553779aa46adb4f2010e33f3def847c7d89b51a8462b227605f6c920fd558a6daf64bc98682e508ae960c" -"0c571870e603ba1fce0c13d53176f353fd319959e13db93eae1359f06e3dd4767c04f824cf34ec7bf8f60161ba1a615db8" -"2852eca9e3869afa711ab9a090660b0dc6cfbea310dda77e02310fbaeacd2636f975838c2dbcdbe9ac2cd85cee28f5e3f0" -"c73abf62f9fa02cd79a7606b7ba855db68a07848b057c3aaf38f1a70086e14616f6f88305a1f9ce6b41378a620d4db3e0e" -"7e1d421590dccaeff86212e232eeb5eb8a8d33a8c9b25ae88f3a7bd5032b4efa68f8af3186a02ffcbf5456f12beccace94" -"c81c360cc4a0dcc642b59f991eec68c59af78139ca60b96d6a18e9535f8995e89bd2cf6a0aef3acffd33d1c0c1b79b6641" -"4a91d9f65b2b4ec65844b96f725d2b4b0c309f3eb9d714e9dd939bbdfd85ce8fb43679aeab13f6c29549949503c9466dbd" -"337c4cdde46d6eacd15f21f4d8fdeaa627a47884c88a9c85f0b731d271a8ea7cb9e04a4a149c23c10f56b3a0476dc77a99" -"9d6e4f813e4b0f805e2a693e2ae4ae0ecc423c9ba5d17b42e691abf83784a582f2b1fd85d1e0a27ba38a500963568b2450" -"363d2c5e3f7b8ba3e5b56e4e9f745a3a710bf2ae233c303068c532ce78ff031e6ab28b705dd94d7db4500909edb5626b8c" -"9bd5ff4f0b4741388f0b91563ee516934c013e901572cba005ac5c535f4f107903be9af7b2793dfb61b5070facbe71eefe" -"1b5600f975c8c38c3a2350d78beadfecb78e981164ae8bc866e732972d3ceef4aac68e15861f9b881d9b51b4edece150bc" -"124b07645defb4202ef5d0e0962db98cae6ed459561c93c74c20bd64362e4f4fffc389a6cd80514604ff22eecc10c9cbc7" -"981d19a8102b24146354c463107c9dc070e29e70df3578022acf72289ef071ab9f9402a544d0399f1b1e5f206b6d46d445" -"f6d612a490e72918e00c853eda8493bef511149e80c9ab56e8b4b8cba3987249f77d060e61760e5792ac321c987c03c260" -"6e9393a7970212992cdbd16448078d5039d4c2c3199714f53278f4f7b1d2e514cf95bdfc078b8bb0db659cb2c3f5cc0289" -"0ea84f05d414c88d2db9e9f8455659b9fa6254405317245fa070d6970cafb4dadb2522b490a5c8e02fe973a8cdbfbfbdbf" -"b01535099ffba3d3896bc4d1189fc570c3e6fdc6469265b8da912772e75dd62ab71be507f700d56cac5e68fd6b57ec1661" -"68ab5258a69625c142a5b1b3519f94be1bde5e51d3bd8ea0c12d5af2fe4615b1b7bd4a96628a4fabc65925ff09718f63bb" -"ebaad98f89bd9543a27b3ff3b5d8bfa89f941a5eb8cc005ccd4a705190e1c9dc6a9f4264e5ee658520a4438e92de854bff" -"c39f8dc7dfbb5de4f14ba63ea16a37d14a7b4610f95b6cffd55e4679b29cedbdf20e7bd16da822fad910c359ee3a68e48a" -"ae6e769b0e291d5d3aa3e2ca9d8d23abe8a1d5349f4991e9300852cc0befb20c2fc0d169306b260763344024f8092cbcc2" -"4c6807363e9fc548a30d5faab3a94b2af0782a2942be80c45d8b0587efd587394ef33c33022436e285806ddffdd32fe363" -"45c3c38ed8d680abeb7a028b44ee6f94d060a14c7019bb6af1f1b5f0a562957d19826d8cc216f9b908c989ccd5415e3525" -"dfe9422ffb5b50b7cc3083dc325544751e5683535d7439d3da2b0bb73bea551dd99e04e0e793804f4774eb6b1daf781d9c" -"aa5128274e599e847862fe309027813d3e4eda0bbeb7201856a5c5d8370e44dabff0bb229c723ba0a6bcf29c44536147de" -"11b7835991018100105bd4329217f7386903fe8e7363cd7b3e893244e245e0a187467664c05b0be1fd429722b9b9a5e319" -"8147fad72776e8a63aab9054fa9d259af0198d088d71d132e6068676a8e9ebb0f616b51ee34aac39c2c2221c7112401727" -"0d75ff4a048363c389e04e9b440ad2032a381ac2cfc54f409caa791e65ee4f5d6cd035008f219b88a803a7382ae447bf65" -"a3df2176b25b3b7b67dabe34decd9a1384dc7a003916ca8fbcb29b3ad6fd8eac5bbbaa3bdfa6c6a3ad9427c4f3ed79fea2" -"6e14c8ce5fa3b4f82c5f7b6d2125916753a7b92ce9b46d45";
\ No newline at end of file diff --git a/numpy/random/src/dsfmt/dSFMT-test-gen.c b/numpy/random/src/dsfmt/dSFMT-test-gen.c deleted file mode 100644 index 697a3010a..000000000 --- a/numpy/random/src/dsfmt/dSFMT-test-gen.c +++ /dev/null @@ -1,58 +0,0 @@ -/* - * cl dSFMT-test-gen.c dSFMT.c -DHAVE_SSE2 -DDSFMT_MEXP=19937 /Ox - * - * gcc dSFMT-test-gen.c dSFMT.c -DHAVE_SSE2 -DDSFMT_MEXP=19937 -o dSFMT - */ - -#include <inttypes.h> -#include <stdio.h> - -#include "dSFMT.h" - - -int main(void) { - int i; - double d; - uint64_t *temp; - uint32_t seed = 1UL; - dsfmt_t state; - dsfmt_init_gen_rand(&state, seed); - double out[1000]; - dsfmt_fill_array_close1_open2(&state, out, 1000); - - FILE *fp; - fp = fopen("dSFMT-testset-1.csv", "w"); - if (fp == NULL) { - printf("Couldn't open file\n"); - return -1; - } - fprintf(fp, "seed, %" PRIu32 "\n", seed); - for (i = 0; i < 1000; i++) { - d = out[i]; - temp = (uint64_t *)&d; - fprintf(fp, "%d, %" PRIu64 "\n", i, *temp); - if (i==999) { - printf("%d, %" PRIu64 "\n", i, *temp); - } - } - fclose(fp); - - seed = 123456789UL; - dsfmt_init_gen_rand(&state, seed); - dsfmt_fill_array_close1_open2(&state, out, 1000); - fp = fopen("dSFMT-testset-2.csv", "w"); - if (fp == NULL) { - printf("Couldn't open file\n"); - return -1; - } - fprintf(fp, "seed, %" PRIu32 "\n", seed); - for (i = 0; i < 1000; i++) { - d = out[i]; - temp = (uint64_t *)&d; - fprintf(fp, "%d, %" PRIu64 "\n", i, *temp); - if (i==999) { - printf("%d, %" PRIu64 "\n", i, *temp); - } - } - fclose(fp); -} diff --git a/numpy/random/src/dsfmt/dSFMT.c b/numpy/random/src/dsfmt/dSFMT.c deleted file mode 100644 index 0f122c26c..000000000 --- a/numpy/random/src/dsfmt/dSFMT.c +++ /dev/null @@ -1,626 +0,0 @@ -/** - * @file dSFMT.c - * @brief double precision SIMD-oriented Fast Mersenne Twister (dSFMT) - * based on IEEE 754 format. - * - * @author Mutsuo Saito (Hiroshima University) - * @author Makoto Matsumoto (Hiroshima University) - * - * Copyright (C) 2007,2008 Mutsuo Saito, Makoto Matsumoto and Hiroshima - * University. All rights reserved. - * - * The new BSD License is applied to this software, see LICENSE.txt - */ -#include <stdio.h> -#include <stdlib.h> -#include <string.h> - -#include "dSFMT-params.h" - -#include "dSFMT-common.h" -#include "dSFMT-jump.h" -#include "dSFMT-poly.h" - -#if defined(__cplusplus) -extern "C" { -#endif - -/** dsfmt internal state vector */ -dsfmt_t dsfmt_global_data; -/** dsfmt mexp for check */ -static const int dsfmt_mexp = DSFMT_MEXP; - -/*---------------- - STATIC FUNCTIONS - ----------------*/ -inline static uint32_t ini_func1(uint32_t x); -inline static uint32_t ini_func2(uint32_t x); -inline static void gen_rand_array_c1o2(dsfmt_t *dsfmt, w128_t *array, int size); -inline static void gen_rand_array_c0o1(dsfmt_t *dsfmt, w128_t *array, int size); -inline static void gen_rand_array_o0c1(dsfmt_t *dsfmt, w128_t *array, int size); -inline static void gen_rand_array_o0o1(dsfmt_t *dsfmt, w128_t *array, int size); -inline static int idxof(int i); -static void initial_mask(dsfmt_t *dsfmt); -static void period_certification(dsfmt_t *dsfmt); - -#if defined(HAVE_SSE2) -/** 1 in 64bit for sse2 */ -static const union X128I_T sse2_int_one = {{1, 1}}; -/** 2.0 double for sse2 */ -static const union X128D_T sse2_double_two = {{2.0, 2.0}}; -/** -1.0 double for sse2 */ -static const union X128D_T sse2_double_m_one = {{-1.0, -1.0}}; -#endif - -/** - * This function simulate a 32-bit array index overlapped to 64-bit - * array of LITTLE ENDIAN in BIG ENDIAN machine. - */ -#if defined(DSFMT_BIG_ENDIAN) -inline static int idxof(int i) { return i ^ 1; } -#else -inline static int idxof(int i) { return i; } -#endif - -#if defined(HAVE_SSE2) -/** - * This function converts the double precision floating point numbers which - * distribute uniformly in the range [1, 2) to those which distribute uniformly - * in the range [0, 1). - * @param w 128bit stracture of double precision floating point numbers (I/O) - */ -inline static void convert_c0o1(w128_t *w) { - w->sd = _mm_add_pd(w->sd, sse2_double_m_one.d128); -} - -/** - * This function converts the double precision floating point numbers which - * distribute uniformly in the range [1, 2) to those which distribute uniformly - * in the range (0, 1]. - * @param w 128bit stracture of double precision floating point numbers (I/O) - */ -inline static void convert_o0c1(w128_t *w) { - w->sd = _mm_sub_pd(sse2_double_two.d128, w->sd); -} - -/** - * This function converts the double precision floating point numbers which - * distribute uniformly in the range [1, 2) to those which distribute uniformly - * in the range (0, 1). - * @param w 128bit stracture of double precision floating point numbers (I/O) - */ -inline static void convert_o0o1(w128_t *w) { - w->si = _mm_or_si128(w->si, sse2_int_one.i128); - w->sd = _mm_add_pd(w->sd, sse2_double_m_one.d128); -} -#else /* standard C and altivec */ -/** - * This function converts the double precision floating point numbers which - * distribute uniformly in the range [1, 2) to those which distribute uniformly - * in the range [0, 1). - * @param w 128bit stracture of double precision floating point numbers (I/O) - */ -inline static void convert_c0o1(w128_t *w) { - w->d[0] -= 1.0; - w->d[1] -= 1.0; -} - -/** - * This function converts the double precision floating point numbers which - * distribute uniformly in the range [1, 2) to those which distribute uniformly - * in the range (0, 1]. - * @param w 128bit stracture of double precision floating point numbers (I/O) - */ -inline static void convert_o0c1(w128_t *w) { - w->d[0] = 2.0 - w->d[0]; - w->d[1] = 2.0 - w->d[1]; -} - -/** - * This function converts the double precision floating point numbers which - * distribute uniformly in the range [1, 2) to those which distribute uniformly - * in the range (0, 1). - * @param w 128bit stracture of double precision floating point numbers (I/O) - */ -inline static void convert_o0o1(w128_t *w) { - w->u[0] |= 1; - w->u[1] |= 1; - w->d[0] -= 1.0; - w->d[1] -= 1.0; -} -#endif - -/** - * This function fills the user-specified array with double precision - * floating point pseudorandom numbers of the IEEE 754 format. - * @param dsfmt dsfmt state vector. - * @param array an 128-bit array to be filled by pseudorandom numbers. - * @param size number of 128-bit pseudorandom numbers to be generated. - */ -inline static void gen_rand_array_c1o2(dsfmt_t *dsfmt, w128_t *array, - int size) { - int i, j; - w128_t lung; - - lung = dsfmt->status[DSFMT_N]; - do_recursion(&array[0], &dsfmt->status[0], &dsfmt->status[DSFMT_POS1], &lung); - for (i = 1; i < DSFMT_N - DSFMT_POS1; i++) { - do_recursion(&array[i], &dsfmt->status[i], &dsfmt->status[i + DSFMT_POS1], - &lung); - } - for (; i < DSFMT_N; i++) { - do_recursion(&array[i], &dsfmt->status[i], &array[i + DSFMT_POS1 - DSFMT_N], - &lung); - } - for (; i < size - DSFMT_N; i++) { - do_recursion(&array[i], &array[i - DSFMT_N], - &array[i + DSFMT_POS1 - DSFMT_N], &lung); - } - for (j = 0; j < 2 * DSFMT_N - size; j++) { - dsfmt->status[j] = array[j + size - DSFMT_N]; - } - for (; i < size; i++, j++) { - do_recursion(&array[i], &array[i - DSFMT_N], - &array[i + DSFMT_POS1 - DSFMT_N], &lung); - dsfmt->status[j] = array[i]; - } - dsfmt->status[DSFMT_N] = lung; -} - -/** - * This function fills the user-specified array with double precision - * floating point pseudorandom numbers of the IEEE 754 format. - * @param dsfmt dsfmt state vector. - * @param array an 128-bit array to be filled by pseudorandom numbers. - * @param size number of 128-bit pseudorandom numbers to be generated. - */ -inline static void gen_rand_array_c0o1(dsfmt_t *dsfmt, w128_t *array, - int size) { - int i, j; - w128_t lung; - - lung = dsfmt->status[DSFMT_N]; - do_recursion(&array[0], &dsfmt->status[0], &dsfmt->status[DSFMT_POS1], &lung); - for (i = 1; i < DSFMT_N - DSFMT_POS1; i++) { - do_recursion(&array[i], &dsfmt->status[i], &dsfmt->status[i + DSFMT_POS1], - &lung); - } - for (; i < DSFMT_N; i++) { - do_recursion(&array[i], &dsfmt->status[i], &array[i + DSFMT_POS1 - DSFMT_N], - &lung); - } - for (; i < size - DSFMT_N; i++) { - do_recursion(&array[i], &array[i - DSFMT_N], - &array[i + DSFMT_POS1 - DSFMT_N], &lung); - convert_c0o1(&array[i - DSFMT_N]); - } - for (j = 0; j < 2 * DSFMT_N - size; j++) { - dsfmt->status[j] = array[j + size - DSFMT_N]; - } - for (; i < size; i++, j++) { - do_recursion(&array[i], &array[i - DSFMT_N], - &array[i + DSFMT_POS1 - DSFMT_N], &lung); - dsfmt->status[j] = array[i]; - convert_c0o1(&array[i - DSFMT_N]); - } - for (i = size - DSFMT_N; i < size; i++) { - convert_c0o1(&array[i]); - } - dsfmt->status[DSFMT_N] = lung; -} - -/** - * This function fills the user-specified array with double precision - * floating point pseudorandom numbers of the IEEE 754 format. - * @param dsfmt dsfmt state vector. - * @param array an 128-bit array to be filled by pseudorandom numbers. - * @param size number of 128-bit pseudorandom numbers to be generated. - */ -inline static void gen_rand_array_o0o1(dsfmt_t *dsfmt, w128_t *array, - int size) { - int i, j; - w128_t lung; - - lung = dsfmt->status[DSFMT_N]; - do_recursion(&array[0], &dsfmt->status[0], &dsfmt->status[DSFMT_POS1], &lung); - for (i = 1; i < DSFMT_N - DSFMT_POS1; i++) { - do_recursion(&array[i], &dsfmt->status[i], &dsfmt->status[i + DSFMT_POS1], - &lung); - } - for (; i < DSFMT_N; i++) { - do_recursion(&array[i], &dsfmt->status[i], &array[i + DSFMT_POS1 - DSFMT_N], - &lung); - } - for (; i < size - DSFMT_N; i++) { - do_recursion(&array[i], &array[i - DSFMT_N], - &array[i + DSFMT_POS1 - DSFMT_N], &lung); - convert_o0o1(&array[i - DSFMT_N]); - } - for (j = 0; j < 2 * DSFMT_N - size; j++) { - dsfmt->status[j] = array[j + size - DSFMT_N]; - } - for (; i < size; i++, j++) { - do_recursion(&array[i], &array[i - DSFMT_N], - &array[i + DSFMT_POS1 - DSFMT_N], &lung); - dsfmt->status[j] = array[i]; - convert_o0o1(&array[i - DSFMT_N]); - } - for (i = size - DSFMT_N; i < size; i++) { - convert_o0o1(&array[i]); - } - dsfmt->status[DSFMT_N] = lung; -} - -/** - * This function fills the user-specified array with double precision - * floating point pseudorandom numbers of the IEEE 754 format. - * @param dsfmt dsfmt state vector. - * @param array an 128-bit array to be filled by pseudorandom numbers. - * @param size number of 128-bit pseudorandom numbers to be generated. - */ -inline static void gen_rand_array_o0c1(dsfmt_t *dsfmt, w128_t *array, - int size) { - int i, j; - w128_t lung; - - lung = dsfmt->status[DSFMT_N]; - do_recursion(&array[0], &dsfmt->status[0], &dsfmt->status[DSFMT_POS1], &lung); - for (i = 1; i < DSFMT_N - DSFMT_POS1; i++) { - do_recursion(&array[i], &dsfmt->status[i], &dsfmt->status[i + DSFMT_POS1], - &lung); - } - for (; i < DSFMT_N; i++) { - do_recursion(&array[i], &dsfmt->status[i], &array[i + DSFMT_POS1 - DSFMT_N], - &lung); - } - for (; i < size - DSFMT_N; i++) { - do_recursion(&array[i], &array[i - DSFMT_N], - &array[i + DSFMT_POS1 - DSFMT_N], &lung); - convert_o0c1(&array[i - DSFMT_N]); - } - for (j = 0; j < 2 * DSFMT_N - size; j++) { - dsfmt->status[j] = array[j + size - DSFMT_N]; - } - for (; i < size; i++, j++) { - do_recursion(&array[i], &array[i - DSFMT_N], - &array[i + DSFMT_POS1 - DSFMT_N], &lung); - dsfmt->status[j] = array[i]; - convert_o0c1(&array[i - DSFMT_N]); - } - for (i = size - DSFMT_N; i < size; i++) { - convert_o0c1(&array[i]); - } - dsfmt->status[DSFMT_N] = lung; -} - -/** - * This function represents a function used in the initialization - * by init_by_array - * @param x 32-bit integer - * @return 32-bit integer - */ -static uint32_t ini_func1(uint32_t x) { - return (x ^ (x >> 27)) * (uint32_t)1664525UL; -} - -/** - * This function represents a function used in the initialization - * by init_by_array - * @param x 32-bit integer - * @return 32-bit integer - */ -static uint32_t ini_func2(uint32_t x) { - return (x ^ (x >> 27)) * (uint32_t)1566083941UL; -} - -/** - * This function initializes the internal state array to fit the IEEE - * 754 format. - * @param dsfmt dsfmt state vector. - */ -static void initial_mask(dsfmt_t *dsfmt) { - int i; - uint64_t *psfmt; - - psfmt = &dsfmt->status[0].u[0]; - for (i = 0; i < DSFMT_N * 2; i++) { - psfmt[i] = (psfmt[i] & DSFMT_LOW_MASK) | DSFMT_HIGH_CONST; - } -} - -/** - * This function certificate the period of 2^{SFMT_MEXP}-1. - * @param dsfmt dsfmt state vector. - */ -static void period_certification(dsfmt_t *dsfmt) { - uint64_t pcv[2] = {DSFMT_PCV1, DSFMT_PCV2}; - uint64_t tmp[2]; - uint64_t inner; - int i; -#if (DSFMT_PCV2 & 1) != 1 - int j; - uint64_t work; -#endif - - tmp[0] = (dsfmt->status[DSFMT_N].u[0] ^ DSFMT_FIX1); - tmp[1] = (dsfmt->status[DSFMT_N].u[1] ^ DSFMT_FIX2); - - inner = tmp[0] & pcv[0]; - inner ^= tmp[1] & pcv[1]; - for (i = 32; i > 0; i >>= 1) { - inner ^= inner >> i; - } - inner &= 1; - /* check OK */ - if (inner == 1) { - return; - } - /* check NG, and modification */ -#if (DSFMT_PCV2 & 1) == 1 - dsfmt->status[DSFMT_N].u[1] ^= 1; -#else - for (i = 1; i >= 0; i--) { - work = 1; - for (j = 0; j < 64; j++) { - if ((work & pcv[i]) != 0) { - dsfmt->status[DSFMT_N].u[i] ^= work; - return; - } - work = work << 1; - } - } -#endif - return; -} - -/*---------------- - PUBLIC FUNCTIONS - ----------------*/ -/** - * This function returns the identification string. The string shows - * the Mersenne exponent, and all parameters of this generator. - * @return id string. - */ -const char *dsfmt_get_idstring(void) { return DSFMT_IDSTR; } - -/** - * This function returns the minimum size of array used for \b - * fill_array functions. - * @return minimum size of array used for fill_array functions. - */ -int dsfmt_get_min_array_size(void) { return DSFMT_N64; } - -/** - * This function fills the internal state array with double precision - * floating point pseudorandom numbers of the IEEE 754 format. - * @param dsfmt dsfmt state vector. - */ -void dsfmt_gen_rand_all(dsfmt_t *dsfmt) { - int i; - w128_t lung; - - lung = dsfmt->status[DSFMT_N]; - do_recursion(&dsfmt->status[0], &dsfmt->status[0], &dsfmt->status[DSFMT_POS1], - &lung); - for (i = 1; i < DSFMT_N - DSFMT_POS1; i++) { - do_recursion(&dsfmt->status[i], &dsfmt->status[i], - &dsfmt->status[i + DSFMT_POS1], &lung); - } - for (; i < DSFMT_N; i++) { - do_recursion(&dsfmt->status[i], &dsfmt->status[i], - &dsfmt->status[i + DSFMT_POS1 - DSFMT_N], &lung); - } - dsfmt->status[DSFMT_N] = lung; -} - -/** - * This function generates double precision floating point - * pseudorandom numbers which distribute in the range [1, 2) to the - * specified array[] by one call. The number of pseudorandom numbers - * is specified by the argument \b size, which must be at least (SFMT_MEXP - * / 128) * 2 and a multiple of two. The function - * get_min_array_size() returns this minimum size. The generation by - * this function is much faster than the following fill_array_xxx functions. - * - * For initialization, init_gen_rand() or init_by_array() must be called - * before the first call of this function. This function can not be - * used after calling genrand_xxx functions, without initialization. - * - * @param dsfmt dsfmt state vector. - * @param array an array where pseudorandom numbers are filled - * by this function. The pointer to the array must be "aligned" - * (namely, must be a multiple of 16) in the SIMD version, since it - * refers to the address of a 128-bit integer. In the standard C - * version, the pointer is arbitrary. - * - * @param size the number of 64-bit pseudorandom integers to be - * generated. size must be a multiple of 2, and greater than or equal - * to (SFMT_MEXP / 128) * 2. - * - * @note \b memalign or \b posix_memalign is available to get aligned - * memory. Mac OSX doesn't have these functions, but \b malloc of OSX - * returns the pointer to the aligned memory block. - */ -void dsfmt_fill_array_close1_open2(dsfmt_t *dsfmt, double array[], int size) { - assert(size % 2 == 0); - assert(size >= DSFMT_N64); - gen_rand_array_c1o2(dsfmt, (w128_t *)array, size / 2); -} - -/** - * This function generates double precision floating point - * pseudorandom numbers which distribute in the range (0, 1] to the - * specified array[] by one call. This function is the same as - * fill_array_close1_open2() except the distribution range. - * - * @param dsfmt dsfmt state vector. - * @param array an array where pseudorandom numbers are filled - * by this function. - * @param size the number of pseudorandom numbers to be generated. - * see also \sa fill_array_close1_open2() - */ -void dsfmt_fill_array_open_close(dsfmt_t *dsfmt, double array[], int size) { - assert(size % 2 == 0); - assert(size >= DSFMT_N64); - gen_rand_array_o0c1(dsfmt, (w128_t *)array, size / 2); -} - -/** - * This function generates double precision floating point - * pseudorandom numbers which distribute in the range [0, 1) to the - * specified array[] by one call. This function is the same as - * fill_array_close1_open2() except the distribution range. - * - * @param array an array where pseudorandom numbers are filled - * by this function. - * @param dsfmt dsfmt state vector. - * @param size the number of pseudorandom numbers to be generated. - * see also \sa fill_array_close1_open2() - */ -void dsfmt_fill_array_close_open(dsfmt_t *dsfmt, double array[], int size) { - assert(size % 2 == 0); - assert(size >= DSFMT_N64); - gen_rand_array_c0o1(dsfmt, (w128_t *)array, size / 2); -} - -/** - * This function generates double precision floating point - * pseudorandom numbers which distribute in the range (0, 1) to the - * specified array[] by one call. This function is the same as - * fill_array_close1_open2() except the distribution range. - * - * @param dsfmt dsfmt state vector. - * @param array an array where pseudorandom numbers are filled - * by this function. - * @param size the number of pseudorandom numbers to be generated. - * see also \sa fill_array_close1_open2() - */ -void dsfmt_fill_array_open_open(dsfmt_t *dsfmt, double array[], int size) { - assert(size % 2 == 0); - assert(size >= DSFMT_N64); - gen_rand_array_o0o1(dsfmt, (w128_t *)array, size / 2); -} - -#if defined(__INTEL_COMPILER) -#pragma warning(disable : 981) -#endif -/** - * This function initializes the internal state array with a 32-bit - * integer seed. - * @param dsfmt dsfmt state vector. - * @param seed a 32-bit integer used as the seed. - * @param mexp caller's mersenne expornent - */ -void dsfmt_chk_init_gen_rand(dsfmt_t *dsfmt, uint32_t seed, int mexp) { - int i; - uint32_t *psfmt; - - /* make sure caller program is compiled with the same MEXP */ - if (mexp != dsfmt_mexp) { - fprintf(stderr, "DSFMT_MEXP doesn't match with dSFMT.c\n"); - exit(1); - } - psfmt = &dsfmt->status[0].u32[0]; - psfmt[idxof(0)] = seed; - for (i = 1; i < (DSFMT_N + 1) * 4; i++) { - psfmt[idxof(i)] = - 1812433253UL * (psfmt[idxof(i - 1)] ^ (psfmt[idxof(i - 1)] >> 30)) + i; - } - initial_mask(dsfmt); - period_certification(dsfmt); - dsfmt->idx = DSFMT_N64; -} - -/** - * This function initializes the internal state array, - * with an array of 32-bit integers used as the seeds - * @param dsfmt dsfmt state vector. - * @param init_key the array of 32-bit integers, used as a seed. - * @param key_length the length of init_key. - * @param mexp caller's mersenne expornent - */ -void dsfmt_chk_init_by_array(dsfmt_t *dsfmt, uint32_t init_key[], - int key_length, int mexp) { - int i, j, count; - uint32_t r; - uint32_t *psfmt32; - int lag; - int mid; - int size = (DSFMT_N + 1) * 4; /* pulmonary */ - - /* make sure caller program is compiled with the same MEXP */ - if (mexp != dsfmt_mexp) { - fprintf(stderr, "DSFMT_MEXP doesn't match with dSFMT.c\n"); - exit(1); - } - if (size >= 623) { - lag = 11; - } else if (size >= 68) { - lag = 7; - } else if (size >= 39) { - lag = 5; - } else { - lag = 3; - } - mid = (size - lag) / 2; - - psfmt32 = &dsfmt->status[0].u32[0]; - memset(dsfmt->status, 0x8b, sizeof(dsfmt->status)); - if (key_length + 1 > size) { - count = key_length + 1; - } else { - count = size; - } - r = ini_func1(psfmt32[idxof(0)] ^ psfmt32[idxof(mid % size)] ^ - psfmt32[idxof((size - 1) % size)]); - psfmt32[idxof(mid % size)] += r; - r += key_length; - psfmt32[idxof((mid + lag) % size)] += r; - psfmt32[idxof(0)] = r; - count--; - for (i = 1, j = 0; (j < count) && (j < key_length); j++) { - r = ini_func1(psfmt32[idxof(i)] ^ psfmt32[idxof((i + mid) % size)] ^ - psfmt32[idxof((i + size - 1) % size)]); - psfmt32[idxof((i + mid) % size)] += r; - r += init_key[j] + i; - psfmt32[idxof((i + mid + lag) % size)] += r; - psfmt32[idxof(i)] = r; - i = (i + 1) % size; - } - for (; j < count; j++) { - r = ini_func1(psfmt32[idxof(i)] ^ psfmt32[idxof((i + mid) % size)] ^ - psfmt32[idxof((i + size - 1) % size)]); - psfmt32[idxof((i + mid) % size)] += r; - r += i; - psfmt32[idxof((i + mid + lag) % size)] += r; - psfmt32[idxof(i)] = r; - i = (i + 1) % size; - } - for (j = 0; j < size; j++) { - r = ini_func2(psfmt32[idxof(i)] + psfmt32[idxof((i + mid) % size)] + - psfmt32[idxof((i + size - 1) % size)]); - psfmt32[idxof((i + mid) % size)] ^= r; - r -= i; - psfmt32[idxof((i + mid + lag) % size)] ^= r; - psfmt32[idxof(i)] = r; - i = (i + 1) % size; - } - initial_mask(dsfmt); - period_certification(dsfmt); - dsfmt->idx = DSFMT_N64; -} -#if defined(__INTEL_COMPILER) -#pragma warning(default : 981) -#endif - -#if defined(__cplusplus) -} -#endif - -extern inline double dsfmt_next_double(dsfmt_state *state); - -extern inline uint64_t dsfmt_next64(dsfmt_state *state); - -extern inline uint32_t dsfmt_next32(dsfmt_state *state); - -void dsfmt_jump(dsfmt_state *state) { dSFMT_jump(state->state, poly_128); };
\ No newline at end of file diff --git a/numpy/random/src/dsfmt/dSFMT.h b/numpy/random/src/dsfmt/dSFMT.h deleted file mode 100644 index 75ef5746f..000000000 --- a/numpy/random/src/dsfmt/dSFMT.h +++ /dev/null @@ -1,697 +0,0 @@ -#pragma once -/** - * @file dSFMT.h - * - * @brief double precision SIMD oriented Fast Mersenne Twister(dSFMT) - * pseudorandom number generator based on IEEE 754 format. - * - * @author Mutsuo Saito (Hiroshima University) - * @author Makoto Matsumoto (Hiroshima University) - * - * Copyright (C) 2007, 2008 Mutsuo Saito, Makoto Matsumoto and - * Hiroshima University. All rights reserved. - * Copyright (C) 2012 Mutsuo Saito, Makoto Matsumoto, - * Hiroshima University and The University of Tokyo. - * All rights reserved. - * - * The new BSD License is applied to this software. - * see LICENSE.txt - * - * @note We assume that your system has inttypes.h. If your system - * doesn't have inttypes.h, you have to typedef uint32_t and uint64_t, - * and you have to define PRIu64 and PRIx64 in this file as follows: - * @verbatim - typedef unsigned int uint32_t - typedef unsigned long long uint64_t - #define PRIu64 "llu" - #define PRIx64 "llx" -@endverbatim - * uint32_t must be exactly 32-bit unsigned integer type (no more, no - * less), and uint64_t must be exactly 64-bit unsigned integer type. - * PRIu64 and PRIx64 are used for printf function to print 64-bit - * unsigned int and 64-bit unsigned int in hexadecimal format. - */ - -#ifndef DSFMT_H -#define DSFMT_H -#if defined(__cplusplus) -extern "C" { -#endif - -#include <assert.h> -#include <stdio.h> - -/* Use NumPy config to determine if SSE2 is present */ -#include "numpy/npy_common.h" -#if defined NPY_HAVE_SSE2_INTRINSICS -#define HAVE_SSE2 1 -#endif - -#if !defined(DSFMT_MEXP) -#ifdef __GNUC__ -#warning "DSFMT_MEXP is not defined. I assume DSFMT_MEXP is 19937." -#endif -#define DSFMT_MEXP 19937 -#endif -/*----------------- - BASIC DEFINITIONS - -----------------*/ -/* Mersenne Exponent. The period of the sequence - * is a multiple of 2^DSFMT_MEXP-1. - * #define DSFMT_MEXP 19937 */ -/** DSFMT generator has an internal state array of 128-bit integers, - * and N is its size. */ -#define DSFMT_N ((DSFMT_MEXP - 128) / 104 + 1) -/** N32 is the size of internal state array when regarded as an array - * of 32-bit integers.*/ -#define DSFMT_N32 (DSFMT_N * 4) -/** N64 is the size of internal state array when regarded as an array - * of 64-bit integers.*/ -#define DSFMT_N64 (DSFMT_N * 2) - -#if !defined(DSFMT_BIG_ENDIAN) -#if defined(__BYTE_ORDER) && defined(__BIG_ENDIAN) -#if __BYTE_ORDER == __BIG_ENDIAN -#define DSFMT_BIG_ENDIAN 1 -#endif -#elif defined(_BYTE_ORDER) && defined(_BIG_ENDIAN) -#if _BYTE_ORDER == _BIG_ENDIAN -#define DSFMT_BIG_ENDIAN 1 -#endif -#elif defined(__BYTE_ORDER__) && defined(__BIG_ENDIAN__) -#if __BYTE_ORDER__ == __BIG_ENDIAN__ -#define DSFMT_BIG_ENDIAN 1 -#endif -#elif defined(BYTE_ORDER) && defined(BIG_ENDIAN) -#if BYTE_ORDER == BIG_ENDIAN -#define DSFMT_BIG_ENDIAN 1 -#endif -#elif defined(__BIG_ENDIAN) || defined(_BIG_ENDIAN) || \ - defined(__BIG_ENDIAN__) || defined(BIG_ENDIAN) -#define DSFMT_BIG_ENDIAN 1 -#endif -#endif - -#if defined(DSFMT_BIG_ENDIAN) && defined(__amd64) -#undef DSFMT_BIG_ENDIAN -#endif - -#if defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901L) -#include <inttypes.h> -#elif defined(_MSC_VER) || defined(__BORLANDC__) -#if !defined(DSFMT_UINT32_DEFINED) && !defined(SFMT_UINT32_DEFINED) -typedef unsigned int uint32_t; -typedef unsigned __int64 uint64_t; -#ifndef UINT64_C -#define UINT64_C(v) (v##ui64) -#endif -#define DSFMT_UINT32_DEFINED -#if !defined(inline) && !defined(__cplusplus) -#define inline __forceinline -#endif -#endif -#else -#include <inttypes.h> -#if !defined(inline) && !defined(__cplusplus) -#if defined(__GNUC__) -#define inline __forceinline__ -#else -#define inline -#endif -#endif -#endif - -#ifndef PRIu64 -#if defined(_MSC_VER) || defined(__BORLANDC__) -#define PRIu64 "I64u" -#define PRIx64 "I64x" -#else -#define PRIu64 "llu" -#define PRIx64 "llx" -#endif -#endif - -#ifndef UINT64_C -#define UINT64_C(v) (v##ULL) -#endif - -/*------------------------------------------ - 128-bit SIMD like data type for standard C - ------------------------------------------*/ -#if defined(HAVE_ALTIVEC) -#if !defined(__APPLE__) -#include <altivec.h> -#endif -/** 128-bit data structure */ -union W128_T { - vector unsigned int s; - uint64_t u[2]; - uint32_t u32[4]; - double d[2]; -}; - -#elif defined(HAVE_SSE2) -#include <emmintrin.h> - -/** 128-bit data structure */ -union W128_T { - __m128i si; - __m128d sd; - uint64_t u[2]; - uint32_t u32[4]; - double d[2]; -}; -#else /* standard C */ -/** 128-bit data structure */ -union W128_T { - uint64_t u[2]; - uint32_t u32[4]; - double d[2]; -}; -#endif - -/** 128-bit data type */ -typedef union W128_T w128_t; - -/** the 128-bit internal state array */ -struct DSFMT_T { - w128_t status[DSFMT_N + 1]; - int idx; -}; -typedef struct DSFMT_T dsfmt_t; - -/** dsfmt internal state vector */ -extern dsfmt_t dsfmt_global_data; -/** dsfmt mexp for check */ -extern const int dsfmt_global_mexp; - -void dsfmt_gen_rand_all(dsfmt_t *dsfmt); -void dsfmt_fill_array_open_close(dsfmt_t *dsfmt, double array[], int size); -void dsfmt_fill_array_close_open(dsfmt_t *dsfmt, double array[], int size); -void dsfmt_fill_array_open_open(dsfmt_t *dsfmt, double array[], int size); -void dsfmt_fill_array_close1_open2(dsfmt_t *dsfmt, double array[], int size); -void dsfmt_chk_init_gen_rand(dsfmt_t *dsfmt, uint32_t seed, int mexp); -void dsfmt_chk_init_by_array(dsfmt_t *dsfmt, uint32_t init_key[], - int key_length, int mexp); -const char *dsfmt_get_idstring(void); -int dsfmt_get_min_array_size(void); - -#if defined(__GNUC__) -#define DSFMT_PRE_INLINE inline static -#define DSFMT_PST_INLINE __attribute__((always_inline)) -#elif defined(_MSC_VER) && _MSC_VER >= 1200 -#define DSFMT_PRE_INLINE __forceinline static -#define DSFMT_PST_INLINE -#else -#define DSFMT_PRE_INLINE inline static -#define DSFMT_PST_INLINE -#endif -DSFMT_PRE_INLINE uint32_t dsfmt_genrand_uint32(dsfmt_t *dsfmt) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE double -dsfmt_genrand_close1_open2(dsfmt_t *dsfmt) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE double -dsfmt_genrand_close_open(dsfmt_t *dsfmt) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE double -dsfmt_genrand_open_close(dsfmt_t *dsfmt) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE double -dsfmt_genrand_open_open(dsfmt_t *dsfmt) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE uint32_t dsfmt_gv_genrand_uint32(void) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE double dsfmt_gv_genrand_close1_open2(void) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE double dsfmt_gv_genrand_close_open(void) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE double dsfmt_gv_genrand_open_close(void) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE double dsfmt_gv_genrand_open_open(void) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE void dsfmt_gv_fill_array_open_close(double array[], - int size) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE void dsfmt_gv_fill_array_close_open(double array[], - int size) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE void dsfmt_gv_fill_array_open_open(double array[], - int size) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE void -dsfmt_gv_fill_array_close1_open2(double array[], int size) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE void dsfmt_gv_init_gen_rand(uint32_t seed) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE void dsfmt_gv_init_by_array(uint32_t init_key[], - int key_length) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE void dsfmt_init_gen_rand(dsfmt_t *dsfmt, - uint32_t seed) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE void dsfmt_init_by_array(dsfmt_t *dsfmt, uint32_t init_key[], - int key_length) DSFMT_PST_INLINE; - -/** - * This function generates and returns unsigned 32-bit integer. - * This is slower than SFMT, only for convenience usage. - * dsfmt_init_gen_rand() or dsfmt_init_by_array() must be called - * before this function. - * @param dsfmt dsfmt internal state date - * @return double precision floating point pseudorandom number - */ -inline static uint32_t dsfmt_genrand_uint32(dsfmt_t *dsfmt) { - uint32_t r; - uint64_t *psfmt64 = &dsfmt->status[0].u[0]; - - if (dsfmt->idx >= DSFMT_N64) { - dsfmt_gen_rand_all(dsfmt); - dsfmt->idx = 0; - } - r = psfmt64[dsfmt->idx++] & 0xffffffffU; - return r; -} - -/** - * This function generates and returns double precision pseudorandom - * number which distributes uniformly in the range [1, 2). This is - * the primitive and faster than generating numbers in other ranges. - * dsfmt_init_gen_rand() or dsfmt_init_by_array() must be called - * before this function. - * @param dsfmt dsfmt internal state date - * @return double precision floating point pseudorandom number - */ -inline static double dsfmt_genrand_close1_open2(dsfmt_t *dsfmt) { - double r; - double *psfmt64 = &dsfmt->status[0].d[0]; - - if (dsfmt->idx >= DSFMT_N64) { - dsfmt_gen_rand_all(dsfmt); - dsfmt->idx = 0; - } - r = psfmt64[dsfmt->idx++]; - return r; -} - -/** - * This function generates and returns unsigned 32-bit integer. - * This is slower than SFMT, only for convenience usage. - * dsfmt_gv_init_gen_rand() or dsfmt_gv_init_by_array() must be called - * before this function. This function uses \b global variables. - * @return double precision floating point pseudorandom number - */ -inline static uint32_t dsfmt_gv_genrand_uint32(void) { - return dsfmt_genrand_uint32(&dsfmt_global_data); -} - -/** - * This function generates and returns double precision pseudorandom - * number which distributes uniformly in the range [1, 2). - * dsfmt_gv_init_gen_rand() or dsfmt_gv_init_by_array() must be called - * before this function. This function uses \b global variables. - * @return double precision floating point pseudorandom number - */ -inline static double dsfmt_gv_genrand_close1_open2(void) { - return dsfmt_genrand_close1_open2(&dsfmt_global_data); -} - -/** - * This function generates and returns double precision pseudorandom - * number which distributes uniformly in the range [0, 1). - * dsfmt_init_gen_rand() or dsfmt_init_by_array() must be called - * before this function. - * @param dsfmt dsfmt internal state date - * @return double precision floating point pseudorandom number - */ -inline static double dsfmt_genrand_close_open(dsfmt_t *dsfmt) { - return dsfmt_genrand_close1_open2(dsfmt) - 1.0; -} - -/** - * This function generates and returns double precision pseudorandom - * number which distributes uniformly in the range [0, 1). - * dsfmt_gv_init_gen_rand() or dsfmt_gv_init_by_array() must be called - * before this function. This function uses \b global variables. - * @return double precision floating point pseudorandom number - */ -inline static double dsfmt_gv_genrand_close_open(void) { - return dsfmt_gv_genrand_close1_open2() - 1.0; -} - -/** - * This function generates and returns double precision pseudorandom - * number which distributes uniformly in the range (0, 1]. - * dsfmt_init_gen_rand() or dsfmt_init_by_array() must be called - * before this function. - * @param dsfmt dsfmt internal state date - * @return double precision floating point pseudorandom number - */ -inline static double dsfmt_genrand_open_close(dsfmt_t *dsfmt) { - return 2.0 - dsfmt_genrand_close1_open2(dsfmt); -} - -/** - * This function generates and returns double precision pseudorandom - * number which distributes uniformly in the range (0, 1]. - * dsfmt_gv_init_gen_rand() or dsfmt_gv_init_by_array() must be called - * before this function. This function uses \b global variables. - * @return double precision floating point pseudorandom number - */ -inline static double dsfmt_gv_genrand_open_close(void) { - return 2.0 - dsfmt_gv_genrand_close1_open2(); -} - -/** - * This function generates and returns double precision pseudorandom - * number which distributes uniformly in the range (0, 1). - * dsfmt_init_gen_rand() or dsfmt_init_by_array() must be called - * before this function. - * @param dsfmt dsfmt internal state date - * @return double precision floating point pseudorandom number - */ -inline static double dsfmt_genrand_open_open(dsfmt_t *dsfmt) { - double *dsfmt64 = &dsfmt->status[0].d[0]; - union { - double d; - uint64_t u; - } r; - - if (dsfmt->idx >= DSFMT_N64) { - dsfmt_gen_rand_all(dsfmt); - dsfmt->idx = 0; - } - r.d = dsfmt64[dsfmt->idx++]; - r.u |= 1; - return r.d - 1.0; -} - -/** - * This function generates and returns double precision pseudorandom - * number which distributes uniformly in the range (0, 1). - * dsfmt_gv_init_gen_rand() or dsfmt_gv_init_by_array() must be called - * before this function. This function uses \b global variables. - * @return double precision floating point pseudorandom number - */ -inline static double dsfmt_gv_genrand_open_open(void) { - return dsfmt_genrand_open_open(&dsfmt_global_data); -} - -/** - * This function generates double precision floating point - * pseudorandom numbers which distribute in the range [1, 2) to the - * specified array[] by one call. This function is the same as - * dsfmt_fill_array_close1_open2() except that this function uses - * \b global variables. - * @param array an array where pseudorandom numbers are filled - * by this function. - * @param size the number of pseudorandom numbers to be generated. - * see also \sa dsfmt_fill_array_close1_open2() - */ -inline static void dsfmt_gv_fill_array_close1_open2(double array[], int size) { - dsfmt_fill_array_close1_open2(&dsfmt_global_data, array, size); -} - -/** - * This function generates double precision floating point - * pseudorandom numbers which distribute in the range (0, 1] to the - * specified array[] by one call. This function is the same as - * dsfmt_gv_fill_array_close1_open2() except the distribution range. - * This function uses \b global variables. - * @param array an array where pseudorandom numbers are filled - * by this function. - * @param size the number of pseudorandom numbers to be generated. - * see also \sa dsfmt_fill_array_close1_open2() and \sa - * dsfmt_gv_fill_array_close1_open2() - */ -inline static void dsfmt_gv_fill_array_open_close(double array[], int size) { - dsfmt_fill_array_open_close(&dsfmt_global_data, array, size); -} - -/** - * This function generates double precision floating point - * pseudorandom numbers which distribute in the range [0, 1) to the - * specified array[] by one call. This function is the same as - * dsfmt_gv_fill_array_close1_open2() except the distribution range. - * This function uses \b global variables. - * @param array an array where pseudorandom numbers are filled - * by this function. - * @param size the number of pseudorandom numbers to be generated. - * see also \sa dsfmt_fill_array_close1_open2() \sa - * dsfmt_gv_fill_array_close1_open2() - */ -inline static void dsfmt_gv_fill_array_close_open(double array[], int size) { - dsfmt_fill_array_close_open(&dsfmt_global_data, array, size); -} - -/** - * This function generates double precision floating point - * pseudorandom numbers which distribute in the range (0, 1) to the - * specified array[] by one call. This function is the same as - * dsfmt_gv_fill_array_close1_open2() except the distribution range. - * This function uses \b global variables. - * @param array an array where pseudorandom numbers are filled - * by this function. - * @param size the number of pseudorandom numbers to be generated. - * see also \sa dsfmt_fill_array_close1_open2() \sa - * dsfmt_gv_fill_array_close1_open2() - */ -inline static void dsfmt_gv_fill_array_open_open(double array[], int size) { - dsfmt_fill_array_open_open(&dsfmt_global_data, array, size); -} - -/** - * This function initializes the internal state array with a 32-bit - * integer seed. - * @param dsfmt dsfmt state vector. - * @param seed a 32-bit integer used as the seed. - */ -inline static void dsfmt_init_gen_rand(dsfmt_t *dsfmt, uint32_t seed) { - dsfmt_chk_init_gen_rand(dsfmt, seed, DSFMT_MEXP); -} - -/** - * This function initializes the internal state array with a 32-bit - * integer seed. This function uses \b global variables. - * @param seed a 32-bit integer used as the seed. - * see also \sa dsfmt_init_gen_rand() - */ -inline static void dsfmt_gv_init_gen_rand(uint32_t seed) { - dsfmt_init_gen_rand(&dsfmt_global_data, seed); -} - -/** - * This function initializes the internal state array, - * with an array of 32-bit integers used as the seeds. - * @param dsfmt dsfmt state vector - * @param init_key the array of 32-bit integers, used as a seed. - * @param key_length the length of init_key. - */ -inline static void dsfmt_init_by_array(dsfmt_t *dsfmt, uint32_t init_key[], - int key_length) { - dsfmt_chk_init_by_array(dsfmt, init_key, key_length, DSFMT_MEXP); -} - -/** - * This function initializes the internal state array, - * with an array of 32-bit integers used as the seeds. - * This function uses \b global variables. - * @param init_key the array of 32-bit integers, used as a seed. - * @param key_length the length of init_key. - * see also \sa dsfmt_init_by_array() - */ -inline static void dsfmt_gv_init_by_array(uint32_t init_key[], int key_length) { - dsfmt_init_by_array(&dsfmt_global_data, init_key, key_length); -} - -#if !defined(DSFMT_DO_NOT_USE_OLD_NAMES) -DSFMT_PRE_INLINE const char *get_idstring(void) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE int get_min_array_size(void) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE void init_gen_rand(uint32_t seed) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE void init_by_array(uint32_t init_key[], - int key_length) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE double genrand_close1_open2(void) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE double genrand_close_open(void) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE double genrand_open_close(void) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE double genrand_open_open(void) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE void fill_array_open_close(double array[], - int size) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE void fill_array_close_open(double array[], - int size) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE void fill_array_open_open(double array[], - int size) DSFMT_PST_INLINE; -DSFMT_PRE_INLINE void fill_array_close1_open2(double array[], - int size) DSFMT_PST_INLINE; - -/** - * This function is just the same as dsfmt_get_idstring(). - * @return id string. - * see also \sa dsfmt_get_idstring() - */ -inline static const char *get_idstring(void) { return dsfmt_get_idstring(); } - -/** - * This function is just the same as dsfmt_get_min_array_size(). - * @return minimum size of array used for fill_array functions. - * see also \sa dsfmt_get_min_array_size() - */ -inline static int get_min_array_size(void) { - return dsfmt_get_min_array_size(); -} - -/** - * This function is just the same as dsfmt_gv_init_gen_rand(). - * @param seed a 32-bit integer used as the seed. - * see also \sa dsfmt_gv_init_gen_rand(), \sa dsfmt_init_gen_rand(). - */ -inline static void init_gen_rand(uint32_t seed) { - dsfmt_gv_init_gen_rand(seed); -} - -/** - * This function is just the same as dsfmt_gv_init_by_array(). - * @param init_key the array of 32-bit integers, used as a seed. - * @param key_length the length of init_key. - * see also \sa dsfmt_gv_init_by_array(), \sa dsfmt_init_by_array(). - */ -inline static void init_by_array(uint32_t init_key[], int key_length) { - dsfmt_gv_init_by_array(init_key, key_length); -} - -/** - * This function is just the same as dsfmt_gv_genrand_close1_open2(). - * @return double precision floating point number. - * see also \sa dsfmt_genrand_close1_open2() \sa - * dsfmt_gv_genrand_close1_open2() - */ -inline static double genrand_close1_open2(void) { - return dsfmt_gv_genrand_close1_open2(); -} - -/** - * This function is just the same as dsfmt_gv_genrand_close_open(). - * @return double precision floating point number. - * see also \sa dsfmt_genrand_close_open() \sa - * dsfmt_gv_genrand_close_open() - */ -inline static double genrand_close_open(void) { - return dsfmt_gv_genrand_close_open(); -} - -/** - * This function is just the same as dsfmt_gv_genrand_open_close(). - * @return double precision floating point number. - * see also \sa dsfmt_genrand_open_close() \sa - * dsfmt_gv_genrand_open_close() - */ -inline static double genrand_open_close(void) { - return dsfmt_gv_genrand_open_close(); -} - -/** - * This function is just the same as dsfmt_gv_genrand_open_open(). - * @return double precision floating point number. - * see also \sa dsfmt_genrand_open_open() \sa - * dsfmt_gv_genrand_open_open() - */ -inline static double genrand_open_open(void) { - return dsfmt_gv_genrand_open_open(); -} - -/** - * This function is juset the same as dsfmt_gv_fill_array_open_close(). - * @param array an array where pseudorandom numbers are filled - * by this function. - * @param size the number of pseudorandom numbers to be generated. - * see also \sa dsfmt_gv_fill_array_open_close(), \sa - * dsfmt_fill_array_close1_open2(), \sa - * dsfmt_gv_fill_array_close1_open2() - */ -inline static void fill_array_open_close(double array[], int size) { - dsfmt_gv_fill_array_open_close(array, size); -} - -/** - * This function is juset the same as dsfmt_gv_fill_array_close_open(). - * @param array an array where pseudorandom numbers are filled - * by this function. - * @param size the number of pseudorandom numbers to be generated. - * see also \sa dsfmt_gv_fill_array_close_open(), \sa - * dsfmt_fill_array_close1_open2(), \sa - * dsfmt_gv_fill_array_close1_open2() - */ -inline static void fill_array_close_open(double array[], int size) { - dsfmt_gv_fill_array_close_open(array, size); -} - -/** - * This function is juset the same as dsfmt_gv_fill_array_open_open(). - * @param array an array where pseudorandom numbers are filled - * by this function. - * @param size the number of pseudorandom numbers to be generated. - * see also \sa dsfmt_gv_fill_array_open_open(), \sa - * dsfmt_fill_array_close1_open2(), \sa - * dsfmt_gv_fill_array_close1_open2() - */ -inline static void fill_array_open_open(double array[], int size) { - dsfmt_gv_fill_array_open_open(array, size); -} - -/** - * This function is juset the same as dsfmt_gv_fill_array_close1_open2(). - * @param array an array where pseudorandom numbers are filled - * by this function. - * @param size the number of pseudorandom numbers to be generated. - * see also \sa dsfmt_fill_array_close1_open2(), \sa - * dsfmt_gv_fill_array_close1_open2() - */ -inline static void fill_array_close1_open2(double array[], int size) { - dsfmt_gv_fill_array_close1_open2(array, size); -} -#endif /* DSFMT_DO_NOT_USE_OLD_NAMES */ - -#if defined(__cplusplus) -} -#endif - -#endif /* DSFMT_H */ - -union random_val_t { - double d; - uint64_t u64; -}; - -typedef struct s_dsfmt_state { - dsfmt_t *state; - int has_uint32; - uint32_t uinteger; - - double *buffered_uniforms; - int buffer_loc; -} dsfmt_state; - -static inline double dsfmt_next_buffer(dsfmt_state *state) { - if (state->buffer_loc < DSFMT_N64) { - double out = state->buffered_uniforms[state->buffer_loc]; - state->buffer_loc++; - return out; - } - dsfmt_fill_array_close1_open2(state->state, state->buffered_uniforms, - DSFMT_N64); - state->buffer_loc = 1; - return state->buffered_uniforms[0]; -} - -static inline double dsfmt_next_double(dsfmt_state *state) { - return dsfmt_next_buffer(state) - 1.0; -} - -static inline uint64_t dsfmt_next64(dsfmt_state *state) { - /* Discard bottom 16 bits */ - uint64_t out; - union random_val_t rv; - rv.d = dsfmt_next_buffer(state); - out = (rv.u64 >> 16) << 32; - rv.d = dsfmt_next_buffer(state); - out |= (rv.u64 >> 16) & 0xffffffff; - return out; -} - -static inline uint32_t dsfmt_next32(dsfmt_state *state) { - /* Discard bottom 16 bits */ - union random_val_t rv; - rv.d = dsfmt_next_buffer(state); - // uint64_t *out = (uint64_t *)&d; - return (uint32_t)((rv.u64 >> 16) & 0xffffffff); -} - -static inline uint64_t dsfmt_next_raw(dsfmt_state *state) { - union random_val_t rv; - rv.d = dsfmt_next_buffer(state); - return rv.u64; -} - -void dsfmt_jump(dsfmt_state *state);
\ No newline at end of file diff --git a/numpy/random/src/legacy/legacy-distributions.c b/numpy/random/src/legacy/legacy-distributions.c index 99665b370..4741a0352 100644 --- a/numpy/random/src/legacy/legacy-distributions.c +++ b/numpy/random/src/legacy/legacy-distributions.c @@ -1,5 +1,6 @@ #include "legacy-distributions.h" + static NPY_INLINE double legacy_double(aug_bitgen_t *aug_state) { return aug_state->bit_generator->next_double(aug_state->bit_generator->state); } @@ -213,6 +214,113 @@ double legacy_exponential(aug_bitgen_t *aug_state, double scale) { return scale * legacy_standard_exponential(aug_state); } + +static RAND_INT_TYPE random_hypergeometric_hyp(bitgen_t *bitgen_state, + RAND_INT_TYPE good, + RAND_INT_TYPE bad, + RAND_INT_TYPE sample) { + RAND_INT_TYPE d1, k, z; + double d2, u, y; + + d1 = bad + good - sample; + d2 = (double)MIN(bad, good); + + y = d2; + k = sample; + while (y > 0.0) { + u = next_double(bitgen_state); + y -= (RAND_INT_TYPE)floor(u + y / (d1 + k)); + k--; + if (k == 0) + break; + } + z = (RAND_INT_TYPE)(d2 - y); + if (good > bad) + z = sample - z; + return z; +} + +/* D1 = 2*sqrt(2/e) */ +/* D2 = 3 - 2*sqrt(3/e) */ +#define D1 1.7155277699214135 +#define D2 0.8989161620588988 +static RAND_INT_TYPE random_hypergeometric_hrua(bitgen_t *bitgen_state, + RAND_INT_TYPE good, + RAND_INT_TYPE bad, + RAND_INT_TYPE sample) { + RAND_INT_TYPE mingoodbad, maxgoodbad, popsize, m, d9; + double d4, d5, d6, d7, d8, d10, d11; + RAND_INT_TYPE Z; + double T, W, X, Y; + + mingoodbad = MIN(good, bad); + popsize = good + bad; + maxgoodbad = MAX(good, bad); + m = MIN(sample, popsize - sample); + d4 = ((double)mingoodbad) / popsize; + d5 = 1.0 - d4; + d6 = m * d4 + 0.5; + d7 = sqrt((double)(popsize - m) * sample * d4 * d5 / (popsize - 1) + 0.5); + d8 = D1 * d7 + D2; + d9 = (RAND_INT_TYPE)floor((double)(m + 1) * (mingoodbad + 1) / (popsize + 2)); + d10 = (loggam(d9 + 1) + loggam(mingoodbad - d9 + 1) + loggam(m - d9 + 1) + + loggam(maxgoodbad - m + d9 + 1)); + d11 = MIN(MIN(m, mingoodbad) + 1.0, floor(d6 + 16 * d7)); + /* 16 for 16-decimal-digit precision in D1 and D2 */ + + while (1) { + X = next_double(bitgen_state); + Y = next_double(bitgen_state); + W = d6 + d8 * (Y - 0.5) / X; + + /* fast rejection: */ + if ((W < 0.0) || (W >= d11)) + continue; + + Z = (RAND_INT_TYPE)floor(W); + T = d10 - (loggam(Z + 1) + loggam(mingoodbad - Z + 1) + loggam(m - Z + 1) + + loggam(maxgoodbad - m + Z + 1)); + + /* fast acceptance: */ + if ((X * (4.0 - X) - 3.0) <= T) + break; + + /* fast rejection: */ + if (X * (X - T) >= 1) + continue; + /* log(0.0) is ok here, since always accept */ + if (2.0 * log(X) <= T) + break; /* acceptance */ + } + + /* this is a correction to HRUA* by Ivan Frohne in rv.py */ + if (good > bad) + Z = m - Z; + + /* another fix from rv.py to allow sample to exceed popsize/2 */ + if (m < sample) + Z = good - Z; + + return Z; +} +#undef D1 +#undef D2 + +static RAND_INT_TYPE random_hypergeometric_original(bitgen_t *bitgen_state, + RAND_INT_TYPE good, + RAND_INT_TYPE bad, + RAND_INT_TYPE sample) +{ + if (sample > 10) { + return random_hypergeometric_hrua(bitgen_state, good, bad, sample); + } else if (sample > 0) { + return random_hypergeometric_hyp(bitgen_state, good, bad, sample); + } else { + return 0; + } +} + + /* * This is a wrapper function that matches the expected template. In the legacy * generator, all int types are long, so this accepts int64 and then converts @@ -223,12 +331,14 @@ double legacy_exponential(aug_bitgen_t *aug_state, double scale) { */ int64_t legacy_random_hypergeometric(bitgen_t *bitgen_state, int64_t good, int64_t bad, int64_t sample) { - return (int64_t)random_hypergeometric(bitgen_state, (RAND_INT_TYPE)good, - (RAND_INT_TYPE)bad, - (RAND_INT_TYPE)sample); + return (int64_t)random_hypergeometric_original(bitgen_state, + (RAND_INT_TYPE)good, + (RAND_INT_TYPE)bad, + (RAND_INT_TYPE)sample); } - int64_t legacy_random_logseries(bitgen_t *bitgen_state, double p) { + +int64_t legacy_random_logseries(bitgen_t *bitgen_state, double p) { return (int64_t)random_logseries(bitgen_state, p); } diff --git a/numpy/random/src/pcg64/pcg64.h b/numpy/random/src/pcg64/pcg64.h index 20d64f8ab..f8578bfb3 100644 --- a/numpy/random/src/pcg64/pcg64.h +++ b/numpy/random/src/pcg64/pcg64.h @@ -53,6 +53,10 @@ #ifdef _WIN32 #include <stdlib.h> #define inline __forceinline +#if _MSC_VER >= 1900 && _M_AMD64 +#include <intrin.h> +#pragma intrinsic(_umul128) +#endif #endif #if __GNUC_GNU_INLINE__ && !defined(__cplusplus) diff --git a/numpy/random/src/xoshiro256/xoshiro256.c b/numpy/random/src/xoshiro256/xoshiro256.c index f5cda2721..a9ac49aa6 100644 --- a/numpy/random/src/xoshiro256/xoshiro256.c +++ b/numpy/random/src/xoshiro256/xoshiro256.c @@ -6,6 +6,7 @@ worldwide. This software is distributed without any warranty. See <http://creativecommons.org/publicdomain/zero/1.0/>. */ +#include <stddef.h> #include "xoshiro256.h" /* This is xoshiro256** 1.0, our all-purpose, rock-solid generator. It has @@ -29,7 +30,7 @@ extern NPY_INLINE uint32_t xoshiro256_next32(xoshiro256_state *state); void xoshiro256_jump(xoshiro256_state *state) { - int i, b; + size_t i, b; static const uint64_t JUMP[] = {0x180ec6d33cfd0aba, 0xd5a61266f0c9392c, 0xa9582618e03fc9aa, 0x39abdc4529b1661c}; uint64_t s0 = 0; diff --git a/numpy/random/src/xoshiro512/xoshiro512.c b/numpy/random/src/xoshiro512/xoshiro512.c index 9fdbed125..4f4af86f0 100644 --- a/numpy/random/src/xoshiro512/xoshiro512.c +++ b/numpy/random/src/xoshiro512/xoshiro512.c @@ -6,6 +6,7 @@ worldwide. This software is distributed without any warranty. See <http://creativecommons.org/publicdomain/zero/1.0/>. */ +#include <stddef.h> #include "xoshiro512.h" /* This is xoshiro512** 1.0, an all-purpose, rock-solid generator. It has @@ -32,7 +33,7 @@ static uint64_t s_placeholder[8]; void xoshiro512_jump(xoshiro512_state *state) { - int i, b, w; + size_t i, b, w; static const uint64_t JUMP[] = {0x33ed89b6e7a353f9, 0x760083d7955323be, 0x2837f2fbb5f22fae, 0x4b8c5674d309511c, 0xb11ac47a7ba28c25, 0xf1be7667092bcc1c, diff --git a/numpy/random/tests/data/dSFMT-testset-1.csv b/numpy/random/tests/data/dSFMT-testset-1.csv deleted file mode 100644 index 9f3f68aee..000000000 --- a/numpy/random/tests/data/dSFMT-testset-1.csv +++ /dev/null @@ -1,1001 +0,0 @@ -seed, 1 -0, 4607719943344484518 -1, 4611291582431749820 -2, 4609448534071823958 -3, 4611106200950903790 -4, 4609580282971267356 -5, 4609720762470045007 -6, 4607636870422133563 -7, 4611678836133173816 -8, 4610735809705346030 -9, 4608789817543785341 -10, 4610393651369520670 -11, 4611623758339670553 -12, 4608444108616433164 -13, 4607465458608330898 -14, 4610235165931416513 -15, 4608103709635697272 -16, 4608100588970381609 -17, 4608641074646810943 -18, 4608616221842668388 -19, 4607549343270027634 -20, 4610344140656206312 -21, 4610535464314418457 -22, 4611312430982209353 -23, 4607401331162512449 -24, 4609116099984352425 -25, 4611257698187747276 -26, 4607645174989853474 -27, 4608837854111428888 -28, 4608691299413198001 -29, 4608930266190921041 -30, 4608332808182207925 -31, 4608421484748440698 -32, 4608383650014337452 -33, 4611545817730524178 -34, 4611496890204636496 -35, 4610861924086585884 -36, 4608194351754852419 -37, 4610283774852204513 -38, 4609474926233319727 -39, 4608727547276598994 -40, 4609717674488922238 -41, 4609423015229480073 -42, 4609349380679114027 -43, 4610751388265382808 -44, 4607733281112877867 -45, 4610591846531581004 -46, 4610735578621585065 -47, 4611251806424951188 -48, 4609460646876881407 -49, 4611551161354804803 -50, 4611249469634822485 -51, 4608821959551777081 -52, 4610052971574464565 -53, 4607777094091618436 -54, 4610851094220499210 -55, 4607702337070583856 -56, 4611414385520522876 -57, 4611576460026085858 -58, 4611473034943841786 -59, 4608897639026169673 -60, 4607883926313481761 -61, 4608797023994348324 -62, 4610894444453627505 -63, 4608376332448068935 -64, 4609893849888152414 -65, 4610202941918038529 -66, 4611056233415559924 -67, 4610877544761426248 -68, 4610994780372079388 -69, 4609067797141254153 -70, 4610982998126297925 -71, 4610789518257002641 -72, 4610598903890296442 -73, 4608025197953180744 -74, 4610185556158198401 -75, 4608809010429661575 -76, 4608768447244428523 -77, 4609124503281513660 -78, 4611101360372917414 -79, 4611042322702996474 -80, 4610858523716073186 -81, 4610375945852385997 -82, 4608483154301089968 -83, 4607798251963762687 -84, 4609119835265156007 -85, 4609188963843200093 -86, 4607266275697659547 -87, 4610923287424505039 -88, 4609618008081028402 -89, 4608201696292514469 -90, 4611530915051149734 -91, 4607364802882841466 -92, 4611236492431573193 -93, 4610807042893722656 -94, 4610080772943107179 -95, 4608570487813674878 -96, 4610901704116296546 -97, 4610535849683417477 -98, 4609487776344637464 -99, 4607684048114742415 -100, 4609756145565247020 -101, 4609839920313018364 -102, 4607300199263955605 -103, 4608349701001143604 -104, 4608563526536655950 -105, 4608211471894692660 -106, 4611561832400682724 -107, 4607607690303987255 -108, 4609884829632656291 -109, 4610554863807344395 -110, 4609913857729390656 -111, 4608533735726217242 -112, 4611205573734963697 -113, 4608825835576863771 -114, 4611237004326669392 -115, 4609254174450399374 -116, 4607622879864982635 -117, 4610819003887309229 -118, 4607698308712315302 -119, 4607833259972583494 -120, 4607681632288287295 -121, 4607910755187354207 -122, 4611122191561013257 -123, 4611656235846796861 -124, 4610565432604340253 -125, 4610908827351842030 -126, 4611268340320100590 -127, 4609313851656821466 -128, 4608528427564169020 -129, 4608200570336395038 -130, 4610538403125154877 -131, 4608888436053455398 -132, 4611543165544411944 -133, 4610250257231662522 -134, 4609297917742740245 -135, 4610400815978769483 -136, 4610528657337572565 -137, 4610705415667566712 -138, 4611077904559980731 -139, 4609927849369442287 -140, 4609397759777226314 -141, 4609010898016502992 -142, 4609120772081864228 -143, 4607237141609153161 -144, 4610571899308981431 -145, 4609008993363640801 -146, 4610287560785901157 -147, 4608623756174376020 -148, 4608653243870415213 -149, 4609547819331301107 -150, 4609673847833623842 -151, 4608235666485521661 -152, 4610114805542994107 -153, 4610479707694856806 -154, 4609103555741806169 -155, 4608500007602890176 -156, 4610672856146553685 -157, 4611610717675563449 -158, 4610730068759051595 -159, 4608685224971234052 -160, 4609413185190097208 -161, 4608897450787802174 -162, 4607186845514749512 -163, 4607348569795042779 -164, 4611009214373065786 -165, 4609635836233645619 -166, 4610977547720639382 -167, 4610180471225090534 -168, 4608711346241655575 -169, 4608627118356010890 -170, 4611222988628267196 -171, 4610585467803239650 -172, 4608953295299223634 -173, 4610267719343621921 -174, 4608384815012305523 -175, 4610760760005158447 -176, 4609888294496573856 -177, 4608487712307397941 -178, 4610046813921597302 -179, 4607901372889700267 -180, 4611648045743259157 -181, 4611287922111507955 -182, 4609242748471226007 -183, 4607755201380013067 -184, 4611417807283839436 -185, 4609947268752725748 -186, 4610889795243117380 -187, 4610682734427309711 -188, 4609430103755831273 -189, 4610392513282878252 -190, 4608028130931026108 -191, 4609354927353821574 -192, 4607999262827197226 -193, 4608772265877184447 -194, 4611465843018312306 -195, 4607857062574325894 -196, 4610666426888867633 -197, 4610600830920094283 -198, 4609205915408410545 -199, 4607853688411879209 -200, 4607747795491854335 -201, 4608262126937716822 -202, 4611534185487112177 -203, 4610074422869556166 -204, 4609237471248388861 -205, 4609177882315404739 -206, 4607287155741697257 -207, 4611476554139586317 -208, 4609517208333105115 -209, 4607319309983389553 -210, 4608526076922090888 -211, 4610748479257207053 -212, 4609931527304084444 -213, 4608009987636217635 -214, 4610407130291313714 -215, 4608891664598965696 -216, 4610296160139884006 -217, 4611484560328101663 -218, 4611019501503618978 -219, 4610246433435823295 -220, 4608562757718917889 -221, 4607471913719777131 -222, 4611613370049530904 -223, 4608319297993304206 -224, 4607745017639172981 -225, 4609731861577957712 -226, 4609663467149266458 -227, 4609543591294589508 -228, 4608524115958204496 -229, 4611293271186698245 -230, 4610645493162693253 -231, 4607841282620685061 -232, 4608714555579084575 -233, 4608149900642705668 -234, 4609881346718346991 -235, 4609652423006025080 -236, 4610576204477932347 -237, 4608419142720251589 -238, 4609292870725527035 -239, 4607915743890091921 -240, 4610631007573258876 -241, 4611091737638956363 -242, 4610866832498942602 -243, 4608679206316379049 -244, 4611254040795706209 -245, 4608564985000526959 -246, 4609448881222436994 -247, 4611606861023266002 -248, 4608930513670169902 -249, 4607323764955464423 -250, 4607288181905687970 -251, 4610373160370490855 -252, 4608411794278861897 -253, 4610212894218458057 -254, 4610694380708700429 -255, 4609922533346642803 -256, 4609392056590729641 -257, 4609803732845487397 -258, 4608878803170308012 -259, 4611524443328391151 -260, 4608174079771415727 -261, 4607408890033763317 -262, 4607845699541088935 -263, 4611555920103058967 -264, 4607479194613061911 -265, 4607653663534995980 -266, 4610070893479029228 -267, 4611538868870820080 -268, 4608567899678704260 -269, 4608770231079288078 -270, 4610411454914405807 -271, 4607820664883172463 -272, 4610714780751262327 -273, 4607194611952450391 -274, 4609087763682578279 -275, 4608165254026394594 -276, 4609234355157830083 -277, 4609341303623897409 -278, 4607843258265880283 -279, 4609385462693327627 -280, 4610305709185463397 -281, 4607607148427164960 -282, 4608714881122218799 -283, 4609651616092383148 -284, 4609231203213271499 -285, 4611257982347817477 -286, 4610152698091154688 -287, 4608144133423192484 -288, 4610573628352437761 -289, 4608544728261953288 -290, 4610198309467009097 -291, 4610449593868273119 -292, 4610593290392091594 -293, 4609046058809591309 -294, 4611622292530238189 -295, 4610657414068263882 -296, 4611165834719653845 -297, 4610350928332385108 -298, 4611352448760095628 -299, 4609948012409647959 -300, 4610309189747788666 -301, 4607755081207867022 -302, 4610231879852064105 -303, 4607888125265337265 -304, 4609172092280206898 -305, 4608588257784565842 -306, 4607741678424158070 -307, 4609025150930086850 -308, 4609393539064468217 -309, 4610911632660167908 -310, 4610958667007850644 -311, 4611286666277101605 -312, 4609804183682242390 -313, 4611608707704948262 -314, 4609669610113267167 -315, 4607666787328441261 -316, 4607581099645179700 -317, 4610190388918185110 -318, 4610216151360211572 -319, 4608284982213104796 -320, 4609043908200033150 -321, 4610094358636174240 -322, 4607727851642292863 -323, 4607748339680477190 -324, 4610333796608461237 -325, 4611630659133526098 -326, 4611011836822995398 -327, 4611271586336335851 -328, 4608676990846072419 -329, 4610486194528414452 -330, 4607576606876065603 -331, 4607719568700291080 -332, 4608551235063435831 -333, 4611011581816455613 -334, 4608841433006333809 -335, 4609590566174099924 -336, 4610751108417575356 -337, 4609783802139185311 -338, 4610078674097919534 -339, 4608133838483219458 -340, 4609277956691541130 -341, 4608489591310323203 -342, 4608190218485836055 -343, 4611079531841410411 -344, 4608880618309483643 -345, 4608911948674088293 -346, 4611291894381153496 -347, 4608451717289459361 -348, 4608796294882212937 -349, 4608414460486049632 -350, 4607422609721232463 -351, 4608080483385313266 -352, 4607622634488318995 -353, 4609289604313013439 -354, 4609239936097680379 -355, 4608372018123900887 -356, 4610702814261825804 -357, 4611274091783983647 -358, 4611215484124931059 -359, 4608990421919168365 -360, 4609097190054835106 -361, 4610994750415795356 -362, 4611072443902954170 -363, 4608952917718970557 -364, 4608900180943861654 -365, 4608934424322310689 -366, 4609731940535270405 -367, 4610297241526025992 -368, 4608524744985785409 -369, 4610233647335282974 -370, 4609397840502965617 -371, 4609931050226720744 -372, 4607823742164438535 -373, 4607386223202154150 -374, 4611077706407577954 -375, 4608540055157729754 -376, 4610737951147572257 -377, 4610902929285966658 -378, 4611385693211960164 -379, 4607224622379354821 -380, 4609166781986849209 -381, 4608748083335025668 -382, 4610443657454469430 -383, 4610468401119056286 -384, 4610937884247828753 -385, 4609084940396193513 -386, 4611415358123142084 -387, 4610501805353766962 -388, 4607767036448986638 -389, 4607461223165192234 -390, 4608226484980255663 -391, 4611607256659641032 -392, 4609945211367732974 -393, 4609006453263783302 -394, 4610265844375613630 -395, 4607694615392738521 -396, 4608606212547814938 -397, 4610034239111814504 -398, 4610103751968466900 -399, 4611088505838050253 -400, 4608851769231884474 -401, 4610288514235425111 -402, 4608505539036108714 -403, 4611453738759382658 -404, 4611101647329150173 -405, 4607983202842737743 -406, 4607628593913051809 -407, 4610817169808213622 -408, 4610274104936495796 -409, 4607686898188473999 -410, 4611494545938384459 -411, 4609445238317096124 -412, 4609809658023272942 -413, 4610395993443671939 -414, 4609532016275791584 -415, 4610018501692092651 -416, 4608683763430851439 -417, 4608548880896401248 -418, 4610478349829709585 -419, 4607855690965045385 -420, 4609679774972563395 -421, 4609301972366993458 -422, 4609957828433462989 -423, 4611601276026033182 -424, 4610886414042292178 -425, 4610540517589250995 -426, 4609329807459066933 -427, 4611012060649555364 -428, 4611004988464520281 -429, 4610092518739796845 -430, 4608982525313436661 -431, 4610220581774992574 -432, 4608389110412341488 -433, 4610577194017978099 -434, 4607777219986546519 -435, 4608552325706694521 -436, 4609384775042120780 -437, 4610819470183619029 -438, 4608862514454376763 -439, 4608050912492519261 -440, 4609954958938789219 -441, 4611451357502982166 -442, 4607476785936630269 -443, 4611329691194458319 -444, 4610683876885297263 -445, 4608922438780754530 -446, 4607347319284557650 -447, 4610212564213298006 -448, 4607187736152210274 -449, 4607821132969264993 -450, 4610701944842365016 -451, 4609138892484188991 -452, 4607579978932469946 -453, 4608297026731285373 -454, 4609117946354613867 -455, 4609873371753866995 -456, 4609883036162181473 -457, 4610617143057865264 -458, 4609705966715129773 -459, 4609266086686667759 -460, 4611092203109148192 -461, 4607277668644988197 -462, 4610243051742855164 -463, 4611488200475462773 -464, 4610159190694085398 -465, 4607187122077884953 -466, 4609178002227614028 -467, 4607200609172450908 -468, 4607203109970409745 -469, 4610157519627986095 -470, 4608168608616624151 -471, 4607556712879928934 -472, 4610602971628266891 -473, 4607272386871519909 -474, 4609226601189759664 -475, 4608821958178910465 -476, 4610337925540682923 -477, 4607756826141445338 -478, 4610670714123277518 -479, 4609997633318663199 -480, 4610992528318514467 -481, 4610508873379935121 -482, 4610548944839799582 -483, 4608576872646763539 -484, 4611475238517289488 -485, 4609969545809504006 -486, 4611604653736723262 -487, 4610513754499061149 -488, 4610047791400434915 -489, 4610466779122303079 -490, 4609199569907073109 -491, 4611355331378329938 -492, 4609211256613089457 -493, 4611345984656025190 -494, 4609276744577281463 -495, 4610677520254288398 -496, 4611565920468553537 -497, 4608887769347254935 -498, 4607891688277052029 -499, 4611210417809931519 -500, 4609181196197018924 -501, 4608620849445253589 -502, 4610338756450099522 -503, 4610235666137930968 -504, 4610190689620274242 -505, 4608156139098624503 -506, 4610233351376219666 -507, 4611116196066412550 -508, 4611244546095227734 -509, 4608354486449139402 -510, 4608722837522685541 -511, 4607298792335449598 -512, 4610940117180049531 -513, 4609905783847698405 -514, 4611068115688450709 -515, 4609567280627055335 -516, 4609668102454567333 -517, 4608575291283631952 -518, 4608606739858036458 -519, 4609920659405064132 -520, 4609633855399730395 -521, 4607420399082287137 -522, 4607497830797814837 -523, 4608734929795542569 -524, 4611677103224173563 -525, 4609895185651955231 -526, 4608551100268458835 -527, 4608794936863357442 -528, 4608839444940253689 -529, 4609723875823547919 -530, 4609134168731540965 -531, 4610864297289180458 -532, 4609561568240290174 -533, 4609455706988469654 -534, 4610110730269692806 -535, 4607590724900811004 -536, 4609841446856073581 -537, 4607519144944801539 -538, 4610958924924618965 -539, 4608058978781928209 -540, 4608930736822030783 -541, 4610339624224904683 -542, 4611268940884582276 -543, 4611614440252938509 -544, 4610283933065539718 -545, 4610827563929259801 -546, 4610238281320018148 -547, 4609068702417082470 -548, 4609965625349945622 -549, 4610567655464689798 -550, 4609517999871284092 -551, 4608853313183377285 -552, 4608597386123068580 -553, 4608596804275711127 -554, 4608806942254133750 -555, 4611595740982862812 -556, 4610653226348519116 -557, 4610010878229382699 -558, 4611430012536690008 -559, 4608194334909286956 -560, 4609770785529395235 -561, 4609636612234158840 -562, 4610467762650198285 -563, 4611250113292757754 -564, 4611123483515753501 -565, 4610256050464540468 -566, 4611554812085476534 -567, 4609545597507432057 -568, 4610251629953739706 -569, 4608097940038860692 -570, 4608939256004427493 -571, 4609549477949346267 -572, 4607856563525396488 -573, 4608407566119329436 -574, 4610977065049540740 -575, 4608677612836947043 -576, 4611670385382852661 -577, 4609169914628845192 -578, 4608385528780825832 -579, 4608431699759708725 -580, 4610213210579325967 -581, 4607790519129120154 -582, 4611460475578903177 -583, 4611645204412117197 -584, 4611045465835867018 -585, 4610795725227740679 -586, 4607610666986980838 -587, 4607713533366355938 -588, 4608008411737790225 -589, 4607218032541409981 -590, 4610712747455657843 -591, 4607322986186115065 -592, 4608609778168478040 -593, 4609117986895835630 -594, 4608387138944308707 -595, 4609405159006321483 -596, 4609201389487900126 -597, 4610814010656557822 -598, 4610461402205528089 -599, 4608856848982780180 -600, 4610009661369407408 -601, 4609531046728456306 -602, 4608781638378145485 -603, 4611071218907304246 -604, 4607718364365206169 -605, 4610766522014845193 -606, 4610418511682022913 -607, 4611489866910598987 -608, 4611024768525348505 -609, 4608411227740737072 -610, 4608347021514952714 -611, 4607229154687220486 -612, 4609527688395331186 -613, 4608610487126715045 -614, 4610163014754346271 -615, 4610119594096556803 -616, 4609099103543638986 -617, 4607960911715387937 -618, 4610543345562112354 -619, 4611673200269784439 -620, 4607890122556287450 -621, 4610510919142595773 -622, 4611000945873569885 -623, 4609861297670464893 -624, 4607365464269028252 -625, 4610263820456779466 -626, 4608382757430988076 -627, 4608592826360850405 -628, 4607897223655826864 -629, 4608783406301942627 -630, 4610831809342653056 -631, 4610592838071858481 -632, 4607625427481844846 -633, 4610803200293531160 -634, 4607315949328468373 -635, 4609568473332490124 -636, 4608018723588381940 -637, 4610473680670701003 -638, 4608867424437758236 -639, 4607226771732395005 -640, 4607648878101783522 -641, 4608407495823699878 -642, 4609303470297933457 -643, 4607995287639912115 -644, 4610604756706603303 -645, 4608065328364362400 -646, 4607659009213858799 -647, 4609407180393559403 -648, 4610161232799622667 -649, 4608312339248283632 -650, 4611365830215244879 -651, 4609241343071166241 -652, 4607187426157508336 -653, 4611008486844877795 -654, 4609348293209960853 -655, 4611430342690450936 -656, 4610022123422557819 -657, 4610662613803950933 -658, 4610421175429479085 -659, 4609631547889552562 -660, 4609940555785407216 -661, 4609822163096232669 -662, 4608970136612861659 -663, 4609427082274890719 -664, 4608697401879465484 -665, 4611207783165609518 -666, 4611373087590380940 -667, 4610545384528497527 -668, 4607694071454287047 -669, 4607913509258771897 -670, 4607226952976335318 -671, 4611367164497924691 -672, 4610773799850733403 -673, 4608923576905855388 -674, 4610829132227252858 -675, 4611539466506594954 -676, 4607450455252831956 -677, 4607924760556738513 -678, 4609257351177318999 -679, 4607886491020993167 -680, 4607262386448907585 -681, 4608805527475164058 -682, 4608519384875417362 -683, 4609768003609528793 -684, 4607990620996706344 -685, 4608000541499168509 -686, 4607514221391064237 -687, 4610596308708149427 -688, 4608457358343713720 -689, 4611109413177548323 -690, 4609292098957449828 -691, 4608275497070553256 -692, 4609949308659603960 -693, 4610508332440425814 -694, 4610523421224858005 -695, 4611628503654168653 -696, 4608988043865917565 -697, 4609452807254068291 -698, 4611008104380823402 -699, 4609415493514583781 -700, 4608204811849219551 -701, 4608154991732011594 -702, 4609565684575358357 -703, 4607201300980991047 -704, 4611578953897989322 -705, 4608949284388541303 -706, 4608953402339590043 -707, 4611094520261253641 -708, 4611564299263181877 -709, 4611244613212746921 -710, 4607665698546290637 -711, 4609929742865966113 -712, 4608756528459788870 -713, 4608559801324682100 -714, 4611161313083363936 -715, 4610640544822605367 -716, 4610461950314271130 -717, 4608429389531989501 -718, 4610594975443340868 -719, 4610653541215203471 -720, 4610354404384656514 -721, 4611322467270517926 -722, 4609004358268238303 -723, 4610113217342068535 -724, 4607247286313434436 -725, 4607936058322365025 -726, 4607498677954044120 -727, 4607367643972642434 -728, 4610903724213995603 -729, 4608398398619525170 -730, 4609011100867415968 -731, 4609286350498400836 -732, 4610564846286379047 -733, 4610610842418549113 -734, 4609379950548700715 -735, 4608749477629127198 -736, 4609389534628643041 -737, 4609709510894589547 -738, 4609720477301256427 -739, 4610433170873472685 -740, 4607581786915955136 -741, 4610426993537088574 -742, 4609893496842706786 -743, 4608182222083733544 -744, 4607415409292672163 -745, 4608909799371727180 -746, 4609682438519448644 -747, 4609837420608110159 -748, 4607722492204198941 -749, 4608063142644927447 -750, 4611212896211946065 -751, 4610459279330601000 -752, 4610766525803719281 -753, 4610541719260518609 -754, 4608446538192511629 -755, 4608529268885531628 -756, 4607702152237957857 -757, 4608797703031075472 -758, 4607439116134819826 -759, 4608311115301487628 -760, 4611675179452768396 -761, 4608076597967526423 -762, 4611585923502702782 -763, 4611007505903425519 -764, 4610334401882712716 -765, 4611292864862708587 -766, 4610520603991775838 -767, 4610790439348561649 -768, 4608020323209861832 -769, 4609132354146195150 -770, 4611648991029158429 -771, 4608415373761338387 -772, 4611222889759456059 -773, 4610394879407915891 -774, 4611223274533537520 -775, 4611048920373264726 -776, 4611203040226595031 -777, 4608581225592953052 -778, 4607944132899105268 -779, 4610553416357950208 -780, 4609183189134981159 -781, 4610931403284842449 -782, 4609626797792255137 -783, 4608437008274383407 -784, 4608841271194024119 -785, 4609511843950082189 -786, 4608432804080683600 -787, 4607886713946305196 -788, 4610350555892554303 -789, 4611041162152526452 -790, 4608810036185927099 -791, 4609731609025465382 -792, 4608387458587420116 -793, 4608846429123315125 -794, 4610376323596472588 -795, 4609423912646885032 -796, 4609218872949994167 -797, 4611375967003041069 -798, 4609485140993387628 -799, 4607604870717557062 -800, 4609495797464442279 -801, 4611456949409319675 -802, 4610344977769413626 -803, 4610598065942935600 -804, 4608013012863891262 -805, 4610252455143552284 -806, 4607700593028756519 -807, 4610045641566183312 -808, 4609480926180737252 -809, 4610275596338864080 -810, 4607659695464558950 -811, 4607219197645073589 -812, 4608177295501330522 -813, 4611273956331899579 -814, 4610913813005660249 -815, 4608470207120093898 -816, 4608174217124512103 -817, 4610065364926597101 -818, 4607349317207213784 -819, 4607602167222023985 -820, 4607657145979677117 -821, 4611508729708873431 -822, 4607908717595303714 -823, 4609727931398518344 -824, 4609540956592359987 -825, 4610440481396242417 -826, 4611346585992438567 -827, 4611152612229187917 -828, 4610384157247730087 -829, 4610830126611132722 -830, 4610272123470087431 -831, 4607234503905390991 -832, 4610613653079230069 -833, 4609179215008588124 -834, 4608441295871321425 -835, 4608116436734160239 -836, 4607605033373857689 -837, 4610599359267200688 -838, 4611379096030431268 -839, 4609842285031861233 -840, 4611250379332137731 -841, 4608487405142537379 -842, 4607380789043538335 -843, 4609546285413174259 -844, 4608919052624376420 -845, 4611474363794717141 -846, 4611079221421189606 -847, 4607871497110868045 -848, 4608251852430693481 -849, 4611271625089563201 -850, 4608142282722604751 -851, 4610614961087854140 -852, 4611030874745849847 -853, 4609674534508351596 -854, 4610748124279118172 -855, 4610214076525417764 -856, 4608915989016466776 -857, 4611186209375381383 -858, 4609729165373960964 -859, 4608171107224247283 -860, 4608322267844345836 -861, 4611385726702876896 -862, 4607526082606428148 -863, 4609300797912528830 -864, 4607995042613073018 -865, 4609544162095522243 -866, 4607273392907536721 -867, 4610915254133616443 -868, 4608528592480458486 -869, 4611065489354804147 -870, 4610750286707033259 -871, 4609777244768435300 -872, 4610807148417457906 -873, 4607877316209555589 -874, 4610316726265842451 -875, 4608771732565061950 -876, 4611471267762612145 -877, 4607984815868159369 -878, 4608744077489931245 -879, 4611032300367986435 -880, 4609572801223705776 -881, 4607388928679638867 -882, 4610440380910085523 -883, 4611677400759288526 -884, 4608231223120382380 -885, 4609826308636672129 -886, 4610729764513821105 -887, 4608945691565841376 -888, 4608283276108322908 -889, 4611090204591740692 -890, 4610600988861462466 -891, 4608814357404053556 -892, 4611331328900205001 -893, 4610440474296736006 -894, 4607431388306045801 -895, 4610821334292221218 -896, 4608554210663333875 -897, 4609824397495829498 -898, 4607541211343519985 -899, 4608435017263349928 -900, 4607219271691108353 -901, 4608430070452421044 -902, 4609082847439943417 -903, 4610866784520449850 -904, 4608287501071307688 -905, 4609218510235145503 -906, 4608114112360957267 -907, 4609922412275378983 -908, 4608601574612929512 -909, 4608063236537296892 -910, 4610507352144992045 -911, 4610831100303954067 -912, 4610989778846895898 -913, 4611131006465079227 -914, 4607610719457154999 -915, 4610658650342157966 -916, 4607418499615550301 -917, 4610402445180375078 -918, 4610463803051786556 -919, 4610040245423397003 -920, 4610291132556872432 -921, 4610915180727115233 -922, 4607198239226330244 -923, 4610719993407015954 -924, 4608790436210431943 -925, 4611518788065155885 -926, 4609410806758155597 -927, 4610354727542013410 -928, 4609032496183417847 -929, 4607612835462448389 -930, 4609119308314247716 -931, 4610676295665807893 -932, 4610030018059715751 -933, 4610396681520935881 -934, 4610115299841718605 -935, 4610531703556384068 -936, 4607313656834232832 -937, 4607826054856970203 -938, 4609717410497090129 -939, 4609043343821435147 -940, 4607629724646231370 -941, 4611347190635269674 -942, 4607431356324177025 -943, 4609743147159956874 -944, 4608919951624732686 -945, 4608549836830011507 -946, 4609835749585271216 -947, 4610001878208091800 -948, 4607727638454636808 -949, 4608140523490695322 -950, 4610951723878037203 -951, 4609561113218416843 -952, 4607375879120504969 -953, 4610968421496640577 -954, 4608729663137359994 -955, 4611521561048982293 -956, 4607647181466306462 -957, 4608815536941397702 -958, 4611410078681334217 -959, 4610883601143579986 -960, 4609217767853028115 -961, 4610569694955441160 -962, 4608142872889589658 -963, 4609078496262967192 -964, 4610075946790752678 -965, 4607952350453678296 -966, 4610919620741525096 -967, 4611050224420434596 -968, 4608163018441029734 -969, 4611368242936987963 -970, 4607644493316907613 -971, 4611292201819050900 -972, 4610919228494056420 -973, 4607225037781465524 -974, 4609803354294636294 -975, 4610012640039408504 -976, 4610054964621136538 -977, 4609178240405976665 -978, 4607687932449852507 -979, 4609284420963602445 -980, 4609123874172501167 -981, 4608282636137081729 -982, 4609729376153713229 -983, 4611206065074370636 -984, 4609819396180228727 -985, 4610891933717707670 -986, 4608390654319867654 -987, 4610530001352182832 -988, 4608968440000980355 -989, 4611276663454436837 -990, 4609638657758409036 -991, 4610986200094730228 -992, 4610734380577234553 -993, 4609408663096464249 -994, 4609878290485950846 -995, 4607522064640469547 -996, 4610791378999926894 -997, 4607540164715119602 -998, 4609346418539511860 -999, 4611057822391293948 diff --git a/numpy/random/tests/data/dSFMT-testset-2.csv b/numpy/random/tests/data/dSFMT-testset-2.csv deleted file mode 100644 index 2ec2d7a51..000000000 --- a/numpy/random/tests/data/dSFMT-testset-2.csv +++ /dev/null @@ -1,1001 +0,0 @@ -seed, 123456789 -0, 4611184391270596398 -1, 4607748002806405513 -2, 4610950593283603895 -3, 4608090221984835117 -4, 4608620401509188549 -5, 4608933259318787038 -6, 4609013189947794152 -7, 4610981701838014712 -8, 4607749361849391199 -9, 4607937126452632796 -10, 4607787877636103702 -11, 4610534911848167216 -12, 4610758085486021829 -13, 4608376900672555174 -14, 4611475749905869564 -15, 4610499914926478618 -16, 4610900925684457179 -17, 4609184136046668572 -18, 4610480821748015478 -19, 4609898134786890046 -20, 4610017288264709787 -21, 4610945461469807520 -22, 4611377383751342990 -23, 4610945102533068783 -24, 4611669662318553242 -25, 4609191925356202514 -26, 4607369493394647937 -27, 4610171428966243908 -28, 4608117572388953111 -29, 4608266910229155519 -30, 4608559463354581151 -31, 4608291596231703883 -32, 4609993154230249063 -33, 4608158820656759750 -34, 4607825861011031752 -35, 4611091529260321033 -36, 4609803558980700204 -37, 4610892045184692457 -38, 4607844754686911799 -39, 4609239584278564809 -40, 4608975935937490223 -41, 4611186452462884146 -42, 4610644474758878544 -43, 4608188546959932195 -44, 4610317408616093972 -45, 4607827178108179262 -46, 4611275432938764949 -47, 4608511443742655969 -48, 4607200952336491646 -49, 4609727773590041393 -50, 4609137674464229063 -51, 4611529391150713249 -52, 4607446200291229812 -53, 4608618724431091553 -54, 4608231197118968153 -55, 4609848377535763472 -56, 4609769454787363651 -57, 4609093431634418146 -58, 4610941233823434235 -59, 4609301175512188901 -60, 4610884340334734656 -61, 4610063958683836346 -62, 4607897080185028324 -63, 4610305504370096344 -64, 4609499707572786607 -65, 4609874865902334026 -66, 4610351583670684094 -67, 4607292638794148255 -68, 4609412782688385863 -69, 4610390851314330496 -70, 4610340712140767101 -71, 4611109929809487388 -72, 4608431958302562464 -73, 4611491800847619506 -74, 4611664179902108071 -75, 4609243488780021515 -76, 4611114015923144350 -77, 4608524512724098287 -78, 4610450327089558934 -79, 4608720370826377375 -80, 4609572763010883699 -81, 4610381056495781763 -82, 4610005210825690556 -83, 4610341565220388232 -84, 4609034688757121677 -85, 4608823563858030255 -86, 4607971981785212736 -87, 4608714648124969040 -88, 4608968779454607977 -89, 4607272316660152324 -90, 4610748035362895446 -91, 4611319709331049410 -92, 4611390662594804501 -93, 4610054939662414847 -94, 4610394601776085983 -95, 4611424622746948620 -96, 4609395571329163407 -97, 4608527159228147662 -98, 4608379506646233343 -99, 4608996098073467200 -100, 4609633861067625056 -101, 4609387288473331641 -102, 4609369354653726335 -103, 4609305386464755625 -104, 4609085462578529471 -105, 4607213117834450226 -106, 4608676147693465985 -107, 4609894016428056597 -108, 4610808023909042858 -109, 4608152398284126687 -110, 4608864655829389907 -111, 4608853043279974421 -112, 4609998495057236534 -113, 4610971344400740897 -114, 4610947677840469570 -115, 4608888205458648733 -116, 4607306226266885750 -117, 4609476897174960733 -118, 4609298769587075756 -119, 4608046849854182276 -120, 4607709982476583693 -121, 4608905629015127110 -122, 4610027478304152622 -123, 4610378827605636632 -124, 4609168469605184889 -125, 4608424320664524511 -126, 4611224145925927248 -127, 4610598390440508158 -128, 4607548101463088366 -129, 4610445109429344448 -130, 4608982837349247842 -131, 4611526772283530460 -132, 4609215133306156120 -133, 4610747257006687691 -134, 4611305960482481336 -135, 4610602687315818756 -136, 4607839670910202130 -137, 4610527601541091852 -138, 4611462737727028786 -139, 4609212169757271808 -140, 4608447721028771557 -141, 4608224903291976145 -142, 4607298632533980318 -143, 4607975637289743494 -144, 4608776340673956742 -145, 4608578170133208214 -146, 4611177429384019220 -147, 4607665628038835535 -148, 4609531011383000577 -149, 4609360969802432085 -150, 4609908488618542164 -151, 4607366945146869514 -152, 4610136614778041377 -153, 4611168569218164361 -154, 4610462572118833671 -155, 4608070981696376487 -156, 4611058778037833630 -157, 4608159294705821382 -158, 4607530053661949689 -159, 4609136441084496096 -160, 4609331241894336822 -161, 4611313107630037707 -162, 4607957053472625207 -163, 4609675126838719650 -164, 4609482958203648215 -165, 4609691585561697624 -166, 4611475423312731438 -167, 4611412076067906505 -168, 4610429355560523848 -169, 4610625875093760541 -170, 4607241368914269823 -171, 4608589893979475221 -172, 4610073186156188115 -173, 4607291233411155158 -174, 4607675047628278616 -175, 4610869395774400226 -176, 4610232438803745722 -177, 4611554162813131112 -178, 4611642473714833781 -179, 4610563419014907913 -180, 4608459192490850885 -181, 4610780149188594220 -182, 4608960376226045768 -183, 4609069361322693819 -184, 4607256923555182400 -185, 4608474664242579394 -186, 4610207389506744572 -187, 4609112003475072746 -188, 4609240603140788354 -189, 4607525348408117354 -190, 4611149396864007205 -191, 4609114686465130118 -192, 4608849128028589904 -193, 4607248401924217556 -194, 4611387244175010347 -195, 4610567855974092668 -196, 4608169346845917282 -197, 4608490452199719856 -198, 4611203596728963611 -199, 4609287639795895008 -200, 4611614031789088265 -201, 4607236671819565824 -202, 4607466608820858068 -203, 4609639323480231358 -204, 4611254610116912557 -205, 4608807229543638034 -206, 4608381564380368174 -207, 4607428682272410485 -208, 4611229637925604134 -209, 4609236526826003496 -210, 4608889880949328886 -211, 4611135901473627148 -212, 4609923185324027506 -213, 4608879482743843897 -214, 4607662449177713187 -215, 4609980811569949840 -216, 4608275147595190524 -217, 4610799005466054235 -218, 4607667917597769158 -219, 4610185589593486163 -220, 4607529819757965470 -221, 4608839547249506178 -222, 4607706145147011973 -223, 4610988495472186250 -224, 4611534361958731542 -225, 4611641549824093970 -226, 4607316484858856298 -227, 4608641757303921184 -228, 4609375357848069574 -229, 4610565635665894453 -230, 4611147322350665952 -231, 4610071054475069545 -232, 4608886005124134993 -233, 4611384240695070553 -234, 4609556577749744408 -235, 4607688273402525356 -236, 4609395656487625483 -237, 4607920617948366178 -238, 4608233544953726639 -239, 4607736865102992897 -240, 4611554956498550667 -241, 4608735997467283056 -242, 4608499613076219431 -243, 4607926707839352263 -244, 4607349468190181214 -245, 4607855564980078814 -246, 4608566548361033733 -247, 4608689878198670581 -248, 4607485839302113425 -249, 4611493753178685166 -250, 4608566613320387204 -251, 4609743179886038389 -252, 4610508594994820223 -253, 4608995913792958562 -254, 4610248353463386070 -255, 4609788192124211795 -256, 4610619330306161425 -257, 4610873599325465005 -258, 4607324385499645328 -259, 4610611165167596515 -260, 4608006298637673371 -261, 4608540339048264499 -262, 4609631136716349669 -263, 4608685013736282276 -264, 4607363759784022848 -265, 4609492611310929004 -266, 4607780070180818716 -267, 4611531753698196550 -268, 4609227266216837458 -269, 4609211002065400677 -270, 4610668395253295080 -271, 4609134381345731597 -272, 4609382192034225627 -273, 4607208308209034488 -274, 4610579328733327647 -275, 4608921603525338555 -276, 4608290209927931669 -277, 4610866583781415461 -278, 4608182329361248100 -279, 4611648549813945436 -280, 4608601920453704621 -281, 4607406218324299637 -282, 4610748358143595351 -283, 4607437422367397844 -284, 4610299319830347312 -285, 4607992330520188137 -286, 4607658701765777668 -287, 4610721959305012250 -288, 4608971493894533044 -289, 4610010722223631500 -290, 4611050493332836673 -291, 4611164520867402836 -292, 4610619993846650787 -293, 4610600062391983254 -294, 4610986071470687711 -295, 4607815296700712791 -296, 4608678841251990428 -297, 4609887779099788759 -298, 4609503319862290027 -299, 4608809762931362117 -300, 4608037449870401927 -301, 4607755403017924034 -302, 4609087730452781738 -303, 4608773046045154889 -304, 4609803415624001168 -305, 4610998875554212160 -306, 4610380022165388956 -307, 4607984105708776524 -308, 4607847620250154418 -309, 4609666480042052524 -310, 4609307223459772378 -311, 4610669103098622941 -312, 4611513493576426284 -313, 4610110724985187558 -314, 4607584875859460118 -315, 4607466337518526743 -316, 4610953875036984820 -317, 4608473324196281668 -318, 4610528420574205379 -319, 4611218523029715214 -320, 4609404517070225101 -321, 4610679296055932161 -322, 4611602007192673713 -323, 4608768227857799294 -324, 4611351262607349204 -325, 4608656666931918232 -326, 4607814222059811944 -327, 4610377735718844205 -328, 4609693488663627404 -329, 4607234605916181353 -330, 4610438458653690136 -331, 4607691881688829838 -332, 4610084067104393530 -333, 4610193058189981242 -334, 4610500065590109969 -335, 4608182288567589802 -336, 4610884979206264676 -337, 4607934930963198287 -338, 4608198333740812601 -339, 4611615912551444803 -340, 4611091273781746311 -341, 4609878217869378267 -342, 4610329799427547900 -343, 4608946066069950044 -344, 4610517372931467061 -345, 4610173879547218394 -346, 4610768143539619524 -347, 4608251912463490886 -348, 4609138858501301814 -349, 4609537774087558923 -350, 4607501599203475779 -351, 4608820206286486654 -352, 4607594549608867563 -353, 4608928529430502872 -354, 4610326793501581341 -355, 4609216901643916714 -356, 4609921023396761286 -357, 4610188250845345370 -358, 4609056567531193554 -359, 4608042356944953893 -360, 4611153374110275273 -361, 4607652688871602388 -362, 4607736758450185452 -363, 4607772815382776660 -364, 4610793989334300613 -365, 4610810029813744832 -366, 4608713713202824549 -367, 4610555523666319407 -368, 4608933966316349782 -369, 4610847233909664040 -370, 4610569003709254271 -371, 4610141934611190870 -372, 4609800637427386711 -373, 4609531911954538534 -374, 4610018946619778104 -375, 4607563033735657544 -376, 4609466294634090519 -377, 4609110904485970900 -378, 4608802716203741548 -379, 4611231193234792818 -380, 4609853965624850005 -381, 4607407678664700238 -382, 4611560957363283790 -383, 4607258843130776963 -384, 4607438437753792222 -385, 4610880518315386981 -386, 4608724997072138032 -387, 4607896367882266335 -388, 4609466683623316620 -389, 4609649679136642775 -390, 4607572059242669390 -391, 4610690224087953221 -392, 4607212888873300995 -393, 4610115548532567091 -394, 4611204182849533970 -395, 4611480154563209673 -396, 4607313745181304733 -397, 4609677304468142434 -398, 4608230866091821000 -399, 4607916785319391722 -400, 4607735989143160304 -401, 4608364795273033367 -402, 4608202139927885958 -403, 4608897400704372931 -404, 4611267249785141575 -405, 4609988674862878902 -406, 4607825900064550736 -407, 4611018040541037989 -408, 4608438772151688632 -409, 4610422591938237999 -410, 4607217184553988938 -411, 4607633087503746743 -412, 4609394147749351901 -413, 4608101641384193571 -414, 4609733515509206078 -415, 4611489547250433971 -416, 4607834589624331833 -417, 4611349716992792673 -418, 4609707846875238752 -419, 4607311797705362203 -420, 4608945328148355588 -421, 4611273525690510581 -422, 4611458884537996759 -423, 4607997755969685936 -424, 4609048489714323017 -425, 4610334128017869552 -426, 4607485869716832613 -427, 4607547499540098372 -428, 4611447798198333473 -429, 4607207442813565439 -430, 4611108178646490883 -431, 4609758124675924332 -432, 4610269457948568827 -433, 4607360068671694963 -434, 4607781179483110631 -435, 4610840076859630697 -436, 4609605188868326206 -437, 4610833404575495679 -438, 4609202151986229830 -439, 4607653465598307819 -440, 4610341806509732173 -441, 4608937637268370608 -442, 4608846981481205936 -443, 4609890399657918800 -444, 4607475810914216914 -445, 4610779510882657410 -446, 4607200291019787105 -447, 4608763897810030884 -448, 4611030953084521579 -449, 4610608205209840707 -450, 4609901665329352338 -451, 4608229933322773774 -452, 4608306405922059711 -453, 4609402784224466904 -454, 4607797912916831810 -455, 4609320676286567523 -456, 4611203509963612873 -457, 4609443449463211381 -458, 4611201121136708490 -459, 4607891679344035909 -460, 4609295647591940857 -461, 4608699650823090334 -462, 4610113773137160513 -463, 4609644998840868353 -464, 4607236971413190205 -465, 4610986387001985169 -466, 4607686165213831157 -467, 4608006708913412573 -468, 4611617607231087789 -469, 4607950605030537282 -470, 4611312308422726037 -471, 4609920921889730694 -472, 4611272051294701454 -473, 4610732866915233164 -474, 4611475736494024667 -475, 4609129855793761412 -476, 4610896503566695638 -477, 4608983293576256239 -478, 4611337113271775442 -479, 4607264202049306366 -480, 4609273459645222412 -481, 4607686257312802596 -482, 4610552669683473434 -483, 4609573159080816112 -484, 4610109994193793014 -485, 4609104807624348930 -486, 4609056640876615682 -487, 4611233171931551808 -488, 4610700243077601839 -489, 4609689839939656894 -490, 4608154258714850667 -491, 4611519937102265713 -492, 4608524210713510379 -493, 4609408429794931452 -494, 4608727835041307081 -495, 4608363974471195432 -496, 4611053981101408157 -497, 4611244348235020563 -498, 4611215359362792075 -499, 4611323939601701219 -500, 4607339198007393537 -501, 4611192785515763411 -502, 4609520870364372480 -503, 4610305448099707859 -504, 4607627137213702268 -505, 4609512376112901200 -506, 4607188668249670063 -507, 4611507107596430103 -508, 4611290552034620332 -509, 4610948015281142465 -510, 4610082188797301672 -511, 4611154579920165202 -512, 4607910614898084038 -513, 4609111687709912685 -514, 4607756890586988655 -515, 4611478346930052063 -516, 4610271854072480776 -517, 4607666773584055448 -518, 4611269065667018778 -519, 4607229932372594880 -520, 4609361761863029782 -521, 4610810902409829664 -522, 4608310590726885309 -523, 4611549741777094242 -524, 4608905382237807476 -525, 4607539324166606283 -526, 4611302527859497090 -527, 4607673514510851852 -528, 4610239139758062881 -529, 4608296614307074920 -530, 4611131538327332418 -531, 4610491790884660304 -532, 4608012090568842826 -533, 4611145939579689859 -534, 4611569174305843109 -535, 4607548241749347055 -536, 4611302507266314629 -537, 4607334076415859573 -538, 4610759794541675536 -539, 4611562195466283509 -540, 4608064277646826273 -541, 4611362206697199696 -542, 4611267027417975453 -543, 4609817290222129321 -544, 4610075404291128380 -545, 4609555606129743990 -546, 4607220569899493231 -547, 4611584841957177930 -548, 4609037839026191075 -549, 4611594336803497113 -550, 4607225960438616513 -551, 4609362154617705500 -552, 4609887291423254556 -553, 4608541390551696577 -554, 4608696812349818364 -555, 4608371224718817057 -556, 4610715234165102256 -557, 4607906422122850842 -558, 4610831254800690212 -559, 4607810400373332275 -560, 4608705747590604299 -561, 4608938946760670556 -562, 4610310158116436046 -563, 4610355131502528018 -564, 4609768625905121586 -565, 4610143261296345738 -566, 4611431373682787281 -567, 4608146686998001641 -568, 4609198539721817636 -569, 4608916158230506393 -570, 4607654288747635129 -571, 4611682519183492769 -572, 4607197631212679817 -573, 4607299807028695407 -574, 4609116180622479613 -575, 4611019095836572557 -576, 4608581189094026112 -577, 4607488328508280547 -578, 4608587490233232612 -579, 4607245708447615950 -580, 4607189799494915135 -581, 4609348574263949313 -582, 4608021918670812153 -583, 4608172706554967110 -584, 4608811025395016288 -585, 4609364751750743520 -586, 4607844470980185823 -587, 4609405096277516268 -588, 4607748139765213490 -589, 4608512257043070004 -590, 4609962195184017357 -591, 4608461665680660962 -592, 4611127630212845842 -593, 4609686172238940069 -594, 4608777755231651430 -595, 4608284543534209439 -596, 4610868067515254496 -597, 4611535716997037852 -598, 4611319738221220860 -599, 4608658969391651641 -600, 4609452813595548756 -601, 4610236109831493974 -602, 4609938178451088584 -603, 4610331640367617101 -604, 4610901433958649983 -605, 4609766058585980491 -606, 4609222434831315585 -607, 4609778306904942608 -608, 4609448207660443683 -609, 4611299794046339746 -610, 4607801595505703392 -611, 4609594326292439532 -612, 4607668862605395543 -613, 4608245023900457864 -614, 4610578512588843180 -615, 4608185699959219467 -616, 4610904181340375013 -617, 4610647304739305074 -618, 4609795287579987586 -619, 4607960349041110093 -620, 4607703003215776639 -621, 4609403905570407605 -622, 4611233143041131400 -623, 4610530479829073842 -624, 4610679919769197229 -625, 4611448708224350289 -626, 4611445633822299312 -627, 4610496480556319861 -628, 4609555553457224207 -629, 4607626163577357218 -630, 4608595404165123581 -631, 4610510352711119715 -632, 4610203134139830798 -633, 4607550008954478579 -634, 4611603434572420257 -635, 4609780364056746558 -636, 4607295948877799964 -637, 4609867047995237092 -638, 4607936708021896797 -639, 4608897965423418533 -640, 4611287469240086203 -641, 4608515945123070881 -642, 4609851530250371283 -643, 4607577382199018499 -644, 4607744147814966969 -645, 4607260472041943130 -646, 4610683962948666275 -647, 4609625943316701593 -648, 4607251851603159602 -649, 4608016163551470839 -650, 4607202891515091580 -651, 4609099272171658208 -652, 4608510662830783836 -653, 4607744672536335386 -654, 4608142194450948613 -655, 4609476103099505412 -656, 4611399217441119768 -657, 4611495773005281088 -658, 4608815211248586470 -659, 4607337589064315457 -660, 4611394644152964336 -661, 4610812001439064700 -662, 4610702350009793284 -663, 4611075442411386625 -664, 4611077060876180663 -665, 4608164209437610624 -666, 4611368259599962784 -667, 4608333197470863467 -668, 4607183015995911227 -669, 4607199710185468635 -670, 4609413972037912933 -671, 4609234714230829818 -672, 4607739028685645905 -673, 4608232319438231981 -674, 4609333542787352994 -675, 4607657722219109388 -676, 4609193924059916664 -677, 4611141187805060655 -678, 4611068281150742947 -679, 4610549552759132313 -680, 4610085533805630329 -681, 4607232810679281805 -682, 4608493447592041083 -683, 4607355443052807819 -684, 4608410340438808883 -685, 4610315775824782427 -686, 4610312241247357403 -687, 4611287815156776852 -688, 4608076401857758978 -689, 4607457081882300105 -690, 4610908420357480199 -691, 4609797527119137644 -692, 4607351051017728429 -693, 4607618982820305008 -694, 4609846699151054310 -695, 4609389871379854176 -696, 4611243148153910479 -697, 4609270449294231868 -698, 4610832482336321517 -699, 4608101914557495685 -700, 4609128450704503077 -701, 4607351438344234793 -702, 4610010340063776057 -703, 4608461610523881117 -704, 4607869099658377415 -705, 4611211613048598168 -706, 4611196065771110369 -707, 4610515053922650643 -708, 4610096469861694516 -709, 4610477093507778048 -710, 4611547661480689243 -711, 4608438911039690892 -712, 4611311311815318674 -713, 4609279386396407118 -714, 4608222142760880731 -715, 4611613394716251191 -716, 4607603661150022989 -717, 4610135239835120022 -718, 4610929039427992532 -719, 4610757208246529003 -720, 4610920496514785256 -721, 4607326205191641070 -722, 4607938491595237155 -723, 4608585902537439220 -724, 4609326104534891368 -725, 4609325776820376036 -726, 4609693740995539995 -727, 4611329366056096595 -728, 4609303022615335557 -729, 4611512548552170265 -730, 4610404528899365728 -731, 4608023620660481005 -732, 4609431135637339890 -733, 4610767321626117704 -734, 4611106580332635792 -735, 4607433026987401919 -736, 4609580917376189588 -737, 4608816125719706388 -738, 4608380327649573838 -739, 4608700977565012592 -740, 4609148128564608995 -741, 4609631585490496912 -742, 4610745913090661333 -743, 4607498234984630394 -744, 4608367220496728902 -745, 4608365876885021447 -746, 4611537321062599251 -747, 4611238252705917535 -748, 4607525503355262497 -749, 4610601812175940986 -750, 4610145907668011789 -751, 4610384184669464682 -752, 4610699305276533889 -753, 4611440399153628650 -754, 4607963045023571960 -755, 4611498554915678298 -756, 4609015832347581911 -757, 4610795942139040060 -758, 4608894432143218464 -759, 4609704019108678046 -760, 4608168143636007672 -761, 4609566697927636482 -762, 4608690694207868944 -763, 4607746195488024521 -764, 4608350743731006452 -765, 4608442252024570087 -766, 4608428784099249674 -767, 4608941009071857822 -768, 4609298165329524240 -769, 4610447927377989769 -770, 4608304643580688447 -771, 4611265394576506233 -772, 4611210499769545678 -773, 4610114198739241967 -774, 4610653279632780678 -775, 4609515286518383576 -776, 4607984314013723903 -777, 4611541983726033367 -778, 4611393756437132236 -779, 4608968117844197920 -780, 4609367443784351367 -781, 4609488775108334110 -782, 4607529648757616057 -783, 4610676930934349350 -784, 4607750265025461672 -785, 4610373465791644318 -786, 4609305678766837551 -787, 4608947449753189724 -788, 4610366767677719066 -789, 4610439177886004542 -790, 4611242968978180676 -791, 4609370292455902521 -792, 4607754584885122450 -793, 4611224375496789735 -794, 4608921239858925416 -795, 4609513753577022933 -796, 4608075523570985167 -797, 4608608957047081948 -798, 4611273688846153770 -799, 4608394757574873003 -800, 4610377036529664140 -801, 4608600356910393592 -802, 4609667431524003711 -803, 4608601585637259149 -804, 4611533564639785432 -805, 4607510309835958191 -806, 4609651505654903275 -807, 4608166496451053374 -808, 4609515171183335141 -809, 4609776525693204395 -810, 4607696284598399608 -811, 4608607508956363891 -812, 4609695267960623947 -813, 4607576367302408137 -814, 4608741052307396862 -815, 4611095472713646530 -816, 4610161900255157770 -817, 4609145054582502965 -818, 4607410140376051944 -819, 4608126518935915215 -820, 4608269617716261203 -821, 4609477491264110038 -822, 4607463147955504958 -823, 4608999294660391637 -824, 4608694924732427850 -825, 4611156031005634796 -826, 4608453663346634965 -827, 4611380857524502488 -828, 4611362793875369801 -829, 4611632478058955853 -830, 4609434664425350531 -831, 4607606564530411276 -832, 4607391976443208678 -833, 4607762558563019180 -834, 4608554249145639939 -835, 4607806692993216225 -836, 4609510831152869655 -837, 4608164624489904634 -838, 4608455009317767175 -839, 4607280108540066925 -840, 4610080527249430824 -841, 4608840198094196329 -842, 4608916984669714190 -843, 4609771655387294402 -844, 4611351501375292078 -845, 4610356649846014183 -846, 4609861702465798084 -847, 4609335612683847594 -848, 4608963836668425606 -849, 4611448716653608808 -850, 4611618237088472583 -851, 4607650248665393412 -852, 4609477068480641193 -853, 4611408250260317487 -854, 4607799702152927524 -855, 4608984567553844241 -856, 4608966215304179278 -857, 4607599007502108199 -858, 4611197470586031919 -859, 4607738821906038713 -860, 4610174343711771016 -861, 4609411396159113704 -862, 4610528341790372072 -863, 4610621185894682737 -864, 4611164850264296206 -865, 4607500722733965482 -866, 4608747074062289526 -867, 4609587390330409056 -868, 4608013778287410191 -869, 4609438917309909895 -870, 4611359511257377419 -871, 4611161903145694224 -872, 4609908825458581349 -873, 4609974364203149964 -874, 4608056454984693014 -875, 4607485841556578933 -876, 4607689636557505920 -877, 4607225026099434704 -878, 4608918180817633858 -879, 4607389899324828547 -880, 4609528891100730648 -881, 4609347474444270651 -882, 4610604256334495724 -883, 4607717534965049292 -884, 4607416814400338843 -885, 4609568365470566179 -886, 4609490489177847460 -887, 4609959177607409888 -888, 4608249931585238164 -889, 4608374394377617948 -890, 4609359264913370700 -891, 4610789661266619275 -892, 4607881230950036624 -893, 4608163786355022310 -894, 4608830462616805753 -895, 4609531962596587483 -896, 4610555549279318514 -897, 4610008765530009024 -898, 4609509527271380682 -899, 4608445793235798406 -900, 4608895922045956617 -901, 4611496044586314375 -902, 4609855938206283389 -903, 4610584515201059904 -904, 4608185787632733541 -905, 4609925998848840417 -906, 4609746471060930910 -907, 4608322802169846228 -908, 4611668609080045996 -909, 4610918346613262546 -910, 4607487495258046096 -911, 4610091716845110326 -912, 4611060358092721143 -913, 4610617258787020006 -914, 4607968616643301279 -915, 4607216453440634248 -916, 4607683961727519867 -917, 4610192441377241514 -918, 4611340079503032986 -919, 4607737818907905432 -920, 4608040273267030617 -921, 4609075420363483026 -922, 4610025209467938351 -923, 4608669897432477872 -924, 4608611467736828996 -925, 4610963769428151250 -926, 4611230933830803123 -927, 4609892039139108424 -928, 4608322827835753071 -929, 4608048405227745232 -930, 4611336950552458383 -931, 4609990562309176924 -932, 4608539034786829718 -933, 4609715165139430182 -934, 4608805499266985258 -935, 4607728070995330274 -936, 4608780877909747196 -937, 4607569412899178661 -938, 4607268788340312926 -939, 4608510300788384404 -940, 4609202712081615466 -941, 4609583146251705462 -942, 4610981698790205568 -943, 4607925526524476327 -944, 4607793604049723576 -945, 4610915422726587727 -946, 4607690153123448022 -947, 4610957908781080072 -948, 4609688199240625930 -949, 4609195637372175715 -950, 4607455193109906152 -951, 4607614996131060051 -952, 4607821739007708428 -953, 4611432473374206640 -954, 4609331676904204846 -955, 4607810059335115947 -956, 4611077768988065423 -957, 4611510065592294343 -958, 4608753144000455824 -959, 4610618261702230984 -960, 4609478955747078670 -961, 4608250680894683660 -962, 4611056070648131063 -963, 4607756102257795122 -964, 4610370838903190290 -965, 4611412764774525666 -966, 4609100881666906368 -967, 4610119679924928715 -968, 4609686905253473358 -969, 4608711239949443984 -970, 4607839187561408271 -971, 4609413459785445169 -972, 4609209994304368132 -973, 4609118705149046785 -974, 4607291458128247233 -975, 4611161411572838996 -976, 4610256654040673624 -977, 4608882855825268963 -978, 4609049328169514708 -979, 4609651814435298462 -980, 4609304465056789103 -981, 4607682759379096849 -982, 4609946393233090661 -983, 4609946524554590950 -984, 4610880973039636436 -985, 4607217356662986962 -986, 4608230276563898969 -987, 4610664933477117472 -988, 4607562227262100270 -989, 4610133121835039282 -990, 4609071027656845298 -991, 4610444138469204749 -992, 4607185460608050805 -993, 4609895459462574326 -994, 4610016322490782234 -995, 4609380549113996677 -996, 4609371524623560982 -997, 4610108153607631096 -998, 4607489006177078361 -999, 4607632190656691768 diff --git a/numpy/random/tests/test_direct.py b/numpy/random/tests/test_direct.py index 722d5e680..534fe2b6d 100644 --- a/numpy/random/tests/test_direct.py +++ b/numpy/random/tests/test_direct.py @@ -6,7 +6,7 @@ from numpy.testing import (assert_equal, assert_allclose, assert_array_equal, assert_raises) import pytest -from numpy.random import (Generator, MT19937, DSFMT, ThreeFry, PCG32, PCG64, +from numpy.random import (Generator, MT19937, ThreeFry, PCG32, PCG64, Philox, Xoshiro256, Xoshiro512, RandomState) from numpy.random.common import interface @@ -510,84 +510,3 @@ class TestMT19937(Base): bit_generator.state = tup actual = rs.integers(2 ** 16) assert_equal(actual, desired) - - -class TestDSFMT(Base): - @classmethod - def setup_class(cls): - cls.bit_generator = DSFMT - cls.bits = 53 - cls.dtype = np.uint64 - cls.data1 = cls._read_csv(join(pwd, './data/dSFMT-testset-1.csv')) - cls.data2 = cls._read_csv(join(pwd, './data/dSFMT-testset-2.csv')) - cls.seed_error_type = TypeError - cls.invalid_seed_types = [] - cls.invalid_seed_values = [(-1,), np.array([2 ** 33]), - (np.array([2 ** 33, 2 ** 33]),)] - - def test_uniform_double(self): - rs = Generator(self.bit_generator(*self.data1['seed'])) - assert_array_equal(uniform_from_dsfmt(self.data1['data']), - rs.random(1000)) - - rs = Generator(self.bit_generator(*self.data2['seed'])) - assert_equal(uniform_from_dsfmt(self.data2['data']), - rs.random(1000)) - - def test_gauss_inv(self): - n = 25 - rs = RandomState(self.bit_generator(*self.data1['seed'])) - gauss = rs.standard_normal(n) - assert_allclose(gauss, - gauss_from_uint(self.data1['data'], n, 'dsfmt')) - - rs = RandomState(self.bit_generator(*self.data2['seed'])) - gauss = rs.standard_normal(25) - assert_allclose(gauss, - gauss_from_uint(self.data2['data'], n, 'dsfmt')) - - def test_seed_out_of_range_array(self): - # GH #82 - rs = Generator(self.bit_generator(*self.data1['seed'])) - assert_raises(ValueError, rs.bit_generator.seed, - [2 ** (self.bits + 1)]) - assert_raises(ValueError, rs.bit_generator.seed, [-1]) - assert_raises(TypeError, rs.bit_generator.seed, - [2 ** (2 * self.bits + 1)]) - - def test_seed_float(self): - # GH #82 - rs = Generator(self.bit_generator(*self.data1['seed'])) - assert_raises(TypeError, rs.bit_generator.seed, np.pi) - assert_raises(TypeError, rs.bit_generator.seed, -np.pi) - - def test_seed_float_array(self): - # GH #82 - rs = Generator(self.bit_generator(*self.data1['seed'])) - assert_raises(TypeError, rs.bit_generator.seed, np.array([np.pi])) - assert_raises(TypeError, rs.bit_generator.seed, np.array([-np.pi])) - assert_raises(TypeError, rs.bit_generator.seed, - np.array([np.pi, -np.pi])) - assert_raises(TypeError, rs.bit_generator.seed, np.array([0, np.pi])) - assert_raises(TypeError, rs.bit_generator.seed, [np.pi]) - assert_raises(TypeError, rs.bit_generator.seed, [0, np.pi]) - - def test_uniform_float(self): - rs = Generator(self.bit_generator(*self.data1['seed'])) - vals = uniform32_from_uint(self.data1['data'], self.bits) - uniforms = rs.random(len(vals), dtype=np.float32) - assert_allclose(uniforms, vals) - assert_equal(uniforms.dtype, np.float32) - - rs = Generator(self.bit_generator(*self.data2['seed'])) - vals = uniform32_from_uint(self.data2['data'], self.bits) - uniforms = rs.random(len(vals), dtype=np.float32) - assert_allclose(uniforms, vals) - assert_equal(uniforms.dtype, np.float32) - - def test_buffer_reset(self): - rs = Generator(self.bit_generator(*self.data1['seed'])) - rs.random(1) - assert rs.bit_generator.state['buffer_loc'] != 382 - rs.bit_generator.seed(*self.data1['seed']) - assert rs.bit_generator.state['buffer_loc'] == 382 diff --git a/numpy/random/tests/test_generator_mt19937.py b/numpy/random/tests/test_generator_mt19937.py index c6259c42a..44324c9e3 100644 --- a/numpy/random/tests/test_generator_mt19937.py +++ b/numpy/random/tests/test_generator_mt19937.py @@ -897,9 +897,9 @@ class TestRandomDist(object): def test_hypergeometric(self): random.bit_generator.seed(self.seed) actual = random.hypergeometric(10.1, 5.5, 14, size=(3, 2)) - desired = np.array([[10, 10], - [10, 10], - [9, 9]]) + desired = np.array([[9, 9], + [10, 9], + [9, 10]]) assert_array_equal(actual, desired) # Test nbad = 0 @@ -1875,7 +1875,7 @@ class TestBroadcast(object): bad_nsample_one = [-1] bad_nsample_two = [4] hypergeom = random.hypergeometric - desired = np.array([1, 1, 1]) + desired = np.array([0, 0, 1]) self.set_seed() actual = hypergeom(ngood * 3, nbad, nsample) @@ -1906,6 +1906,11 @@ class TestBroadcast(object): assert_raises(ValueError, hypergeom, 10, 10, -1) assert_raises(ValueError, hypergeom, 10, 10, 25) + # ValueError for arguments that are too big. + assert_raises(ValueError, hypergeom, 2**30, 10, 20) + assert_raises(ValueError, hypergeom, 999, 2**31, 50) + assert_raises(ValueError, hypergeom, 999, [2**29, 2**30], 1000) + def test_logseries(self): p = [0.5] bad_p_one = [2] diff --git a/numpy/random/tests/test_generator_mt19937_regressions.py b/numpy/random/tests/test_generator_mt19937_regressions.py index cd4e08d6f..7dca65071 100644 --- a/numpy/random/tests/test_generator_mt19937_regressions.py +++ b/numpy/random/tests/test_generator_mt19937_regressions.py @@ -23,15 +23,8 @@ class TestRegression(object): assert_(np.all(mt19937.hypergeometric(18, 3, 11, size=10) > 0)) # Test for ticket #5623 - args = [ - (2**20 - 2, 2**20 - 2, 2**20 - 2), # Check for 32-bit systems - ] - is_64bits = sys.maxsize > 2**32 - if is_64bits and sys.platform != 'win32': - # Check for 64-bit systems - args.append((2**40 - 2, 2**40 - 2, 2**40 - 2)) - for arg in args: - assert_(mt19937.hypergeometric(*arg) > 0) + args = (2**20 - 2, 2**20 - 2, 2**20 - 2) # Check for 32-bit systems + assert_(mt19937.hypergeometric(*args) > 0) def test_logseries_convergence(self): # Test for ticket #923 diff --git a/numpy/random/tests/test_smoke.py b/numpy/random/tests/test_smoke.py index c94e1a285..ba7f72159 100644 --- a/numpy/random/tests/test_smoke.py +++ b/numpy/random/tests/test_smoke.py @@ -5,7 +5,7 @@ from functools import partial import numpy as np import pytest from numpy.testing import assert_equal, assert_, assert_array_equal -from numpy.random import (Generator, MT19937, DSFMT, ThreeFry, +from numpy.random import (Generator, MT19937, ThreeFry, PCG32, PCG64, Philox, Xoshiro256, Xoshiro512, entropy) @@ -813,18 +813,6 @@ class TestXoshiro512(RNG): cls._extra_setup() -class TestDSFMT(RNG): - @classmethod - def setup_class(cls): - cls.bit_generator = DSFMT - cls.advance = None - cls.seed = [12345] - cls.rg = Generator(cls.bit_generator(*cls.seed)) - cls.initial_state = cls.rg.bit_generator.state - cls._extra_setup() - cls.seed_vector_bits = 32 - - class TestEntropy(object): def test_entropy(self): e1 = entropy.random_entropy() diff --git a/tools/refguide_check.py b/tools/refguide_check.py index 74dbad78b..c20807267 100644 --- a/tools/refguide_check.py +++ b/tools/refguide_check.py @@ -505,7 +505,8 @@ class Checker(doctest.OutputChecker): obj_pattern = re.compile('at 0x[0-9a-fA-F]+>') int_pattern = re.compile('^[0-9]+L?$') vanilla = doctest.OutputChecker() - rndm_markers = {'# random', '# Random', '#random', '#Random', "# may vary"} + rndm_markers = {'# random', '# Random', '#random', '#Random', "# may vary", + "# uninitialized", "#uninitialized"} stopwords = {'plt.', '.hist', '.show', '.ylim', '.subplot(', 'set_title', 'imshow', 'plt.show', '.axis(', '.plot(', '.bar(', '.title', '.ylabel', '.xlabel', 'set_ylim', 'set_xlim', |