diff --git a/CHANGELOG.md b/CHANGELOG.md index 61e54d3..f1d9c31 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,9 +11,15 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Changed * Pinned Cython in the Coverity Scan workflow so generated code stays stable between scans, and added `coverity/README.md` documenting the known Cython-boilerplate false positives and the scan review checklist [gh-164](https://github.com/IntelPython/mkl_random/pull/164) +* Sped up `randint` for power-of-two ranges at or above `INT_MAX` [gh-172](https://github.com/IntelPython/mkl_random/pull/172) +* The random streams for power-of-two `randint` ranges at or above `INT_MAX` have changed: with a fixed seed these now produce different (but equally valid) samples. All other streams are unaffected. [gh-172](https://github.com/IntelPython/mkl_random/pull/172) +* Sped up `randint` for `bool`, `uint8`, `int8`, `uint16` and `int16`; generated values are unchanged [gh-172](https://github.com/IntelPython/mkl_random/pull/172) ### Fixed * Fixed `uniform` to return a Python `float` for scalar bounds with `size=None` instead of a 0-d array [gh-167](https://github.com/IntelPython/mkl_random/pull/167) +* Fixed `randint` returning `high` for `int64`, `uint64` and the default `int` dtype when the range is at or above `INT_MAX` [gh-172](https://github.com/IntelPython/mkl_random/pull/172) +* Fixed the integer fills silently under-filling requests larger than two `MKL_INT_MAX` chunks [gh-172](https://github.com/IntelPython/mkl_random/pull/172) +* Fixed `multinomial` under-filling large outputs by decrementing its chunk counter by elements instead of draws [gh-172](https://github.com/IntelPython/mkl_random/pull/172) ## [1.5.0] (08/12/2026) diff --git a/mkl_random/src/mkl_distributions.cpp b/mkl_random/src/mkl_distributions.cpp index 01cd256..9d0e780 100644 --- a/mkl_random/src/mkl_distributions.cpp +++ b/mkl_random/src/mkl_distributions.cpp @@ -1288,8 +1288,9 @@ void irk_multinomial_vec(irk_state *state, err = viRngMultinomial(VSL_RNG_METHOD_MULTINOMIAL_MULTPOISSON, state->stream, MKL_INT_MAX, res, n, k, pvec); assert(err == VSL_STATUS_OK); + /* len counts draws, res counts ints. */ res += k * MKL_INT_MAX; - len -= k * MKL_INT_MAX; + len -= MKL_INT_MAX; } err = viRngMultinomial(VSL_RNG_METHOD_MULTINOMIAL_MULTPOISSON, @@ -1736,26 +1737,46 @@ void irk_long_vec(irk_state *state, npy_intp len, long *res) res[i] = (long)(ulptr[i] >> 1); } +/* viRngUniform emits 32-bit integers, so a narrower fill needs staging. A + * cache-resident tile avoids allocating 4 bytes per element for the request. */ +template +static void irk_rand_narrow_fill(irk_state *state, + npy_intp len, + T *res, + const int lo, + const int hi) +{ + int err = 0; + const npy_intp tile_len = 4096; + int tile[tile_len]; + + while (len > 0) { + const npy_intp n = (len < tile_len) ? len : tile_len; + npy_intp i = 0; + err = viRngUniform(VSL_RNG_METHOD_UNIFORM_STD, state->stream, (int)n, + tile, lo, hi + 1); + assert(err == VSL_STATUS_OK); + + DIST_PRAGMA_VECTOR + for (i = 0; i < n; ++i) + res[i] = (T)tile[i]; + + res += n; + len -= n; + } +} + void irk_rand_bool_vec(irk_state *state, npy_intp len, npy_bool *res, const npy_bool lo, const npy_bool hi) { - int err = 0; npy_intp i = 0; - int *buf = nullptr; if (len < 1) return; - if (len > MKL_INT_MAX) { - irk_rand_bool_vec(state, MKL_INT_MAX, res, lo, hi); - - res += MKL_INT_MAX; - len -= MKL_INT_MAX; - } - if (lo == hi) { DIST_PRAGMA_VECTOR for (i = 0; i < len; ++i) @@ -1765,18 +1786,7 @@ void irk_rand_bool_vec(irk_state *state, } assert((lo == 0) && (hi == 1)); - buf = (int *)mkl_malloc(len * sizeof(int), 64); - assert(buf != nullptr); - - err = viRngUniform(VSL_RNG_METHOD_UNIFORM_STD, state->stream, len, buf, - (int)lo, (int)hi + 1); - assert(err == VSL_STATUS_OK); - - DIST_PRAGMA_VECTOR - for (i = 0; i < len; ++i) - res[i] = (npy_bool)buf[i]; - - mkl_free(buf); + irk_rand_narrow_fill(state, len, res, (int)lo, (int)hi); } void irk_rand_uint8_vec(irk_state *state, @@ -1785,20 +1795,11 @@ void irk_rand_uint8_vec(irk_state *state, const npy_uint8 lo, const npy_uint8 hi) { - int err = 0; npy_intp i = 0; - int *buf = nullptr; if (len < 1) return; - if (len > MKL_INT_MAX) { - irk_rand_uint8_vec(state, MKL_INT_MAX, res, lo, hi); - - res += MKL_INT_MAX; - len -= MKL_INT_MAX; - } - if (lo == hi) { DIST_PRAGMA_VECTOR for (i = 0; i < len; ++i) @@ -1808,18 +1809,7 @@ void irk_rand_uint8_vec(irk_state *state, } assert(lo < hi); - buf = (int *)mkl_malloc(len * sizeof(int), 64); - assert(buf != nullptr); - - err = viRngUniform(VSL_RNG_METHOD_UNIFORM_STD, state->stream, len, buf, - (int)lo, (int)hi + 1); - assert(err == VSL_STATUS_OK); - - DIST_PRAGMA_VECTOR - for (i = 0; i < len; ++i) - res[i] = (npy_uint8)buf[i]; - - mkl_free(buf); + irk_rand_narrow_fill(state, len, res, (int)lo, (int)hi); } void irk_rand_int8_vec(irk_state *state, @@ -1828,20 +1818,11 @@ void irk_rand_int8_vec(irk_state *state, const npy_int8 lo, const npy_int8 hi) { - int err = 0; npy_intp i = 0; - int *buf = nullptr; if (len < 1) return; - if (len > MKL_INT_MAX) { - irk_rand_int8_vec(state, MKL_INT_MAX, res, lo, hi); - - res += MKL_INT_MAX; - len -= MKL_INT_MAX; - } - if (lo == hi) { DIST_PRAGMA_VECTOR for (i = 0; i < len; ++i) @@ -1851,18 +1832,7 @@ void irk_rand_int8_vec(irk_state *state, } assert(lo < hi); - buf = (int *)mkl_malloc(len * sizeof(int), 64); - assert(buf != nullptr); - - err = viRngUniform(VSL_RNG_METHOD_UNIFORM_STD, state->stream, len, buf, - (int)lo, (int)hi + 1); - assert(err == VSL_STATUS_OK); - - DIST_PRAGMA_VECTOR - for (i = 0; i < len; ++i) - res[i] = (npy_int8)buf[i]; - - mkl_free(buf); + irk_rand_narrow_fill(state, len, res, (int)lo, (int)hi); } void irk_rand_uint16_vec(irk_state *state, @@ -1871,20 +1841,11 @@ void irk_rand_uint16_vec(irk_state *state, const npy_uint16 lo, const npy_uint16 hi) { - int err = 0; npy_intp i = 0; - int *buf = nullptr; if (len < 1) return; - if (len > MKL_INT_MAX) { - irk_rand_uint16_vec(state, MKL_INT_MAX, res, lo, hi); - - res += MKL_INT_MAX; - len -= MKL_INT_MAX; - } - if (lo == hi) { DIST_PRAGMA_VECTOR for (i = 0; i < len; ++i) @@ -1894,18 +1855,7 @@ void irk_rand_uint16_vec(irk_state *state, } assert(lo < hi); - buf = (int *)mkl_malloc(len * sizeof(int), 64); - assert(buf != nullptr); - - err = viRngUniform(VSL_RNG_METHOD_UNIFORM_STD, state->stream, len, buf, - (int)lo, (int)hi + 1); - assert(err == VSL_STATUS_OK); - - DIST_PRAGMA_VECTOR - for (i = 0; i < len; ++i) - res[i] = (npy_uint16)buf[i]; - - mkl_free(buf); + irk_rand_narrow_fill(state, len, res, (int)lo, (int)hi); } void irk_rand_int16_vec(irk_state *state, @@ -1914,20 +1864,11 @@ void irk_rand_int16_vec(irk_state *state, const npy_int16 lo, const npy_int16 hi) { - int err = 0; npy_intp i = 0; - int *buf = nullptr; if (len < 1) return; - if (len > MKL_INT_MAX) { - irk_rand_int16_vec(state, MKL_INT_MAX, res, lo, hi); - - res += MKL_INT_MAX; - len -= MKL_INT_MAX; - } - if (lo == hi) { DIST_PRAGMA_VECTOR for (i = 0; i < len; ++i) @@ -1937,18 +1878,7 @@ void irk_rand_int16_vec(irk_state *state, } assert(lo < hi); - buf = (int *)mkl_malloc(len * sizeof(int), 64); - assert(buf != nullptr); - - err = viRngUniform(VSL_RNG_METHOD_UNIFORM_STD, state->stream, len, buf, - (int)lo, (int)hi + 1); - assert(err == VSL_STATUS_OK); - - DIST_PRAGMA_VECTOR - for (i = 0; i < len; ++i) - res[i] = (npy_int16)buf[i]; - - mkl_free(buf); + irk_rand_narrow_fill(state, len, res, (int)lo, (int)hi); } void irk_rand_uint32_vec(irk_state *state, @@ -1963,7 +1893,7 @@ void irk_rand_uint32_vec(irk_state *state, if (len < 1) return; - if (len > MKL_INT_MAX) { + while (len > MKL_INT_MAX) { irk_rand_uint32_vec(state, MKL_INT_MAX, res, lo, hi); res += MKL_INT_MAX; @@ -2017,7 +1947,7 @@ void irk_rand_int32_vec(irk_state *state, if (len < 1) return; - if (len > MKL_INT_MAX) { + while (len > MKL_INT_MAX) { irk_rand_int32_vec(state, MKL_INT_MAX, res, lo, hi); res += MKL_INT_MAX; @@ -2054,7 +1984,7 @@ void irk_rand_uint64_vec(irk_state *state, if (len < 1) return; - if (len > MKL_INT_MAX) { + while (len > MKL_INT_MAX) { irk_rand_uint64_vec(state, MKL_INT_MAX, res, lo, hi); res += MKL_INT_MAX; @@ -2070,6 +2000,8 @@ void irk_rand_uint64_vec(irk_state *state, return; } + /* Inclusive maximum offset, not the count: the masked branch derives its + * mask from rng, and an incremented rng would admit hi + 1. */ rng = hi - lo; if (!rng) { DIST_PRAGMA_VECTOR @@ -2079,14 +2011,13 @@ void irk_rand_uint64_vec(irk_state *state, return; } - rng++; - - if (rng <= (npy_uint64)INT_MAX) { + if (rng < (npy_uint64)INT_MAX) { int *buf = (int *)mkl_malloc(len * sizeof(int), 64); assert(buf != nullptr); + /* viRngUniform's upper bound is exclusive, hence rng + 1. */ err = viRngUniform(VSL_RNG_METHOD_UNIFORM_STD, state->stream, len, buf, - 0, (int)rng); + 0, (int)(rng + 1)); assert(err == VSL_STATUS_OK); DIST_PRAGMA_VECTOR @@ -2107,26 +2038,58 @@ void irk_rand_uint64_vec(irk_state *state, mask |= mask >> 16; mask |= mask >> 32; - buf = (npy_uint64 *)mkl_malloc(len * sizeof(npy_uint64), 64); - assert(buf != nullptr); + if (mask == rng) { + /* rng + 1 is a power of two, so masking alone confines every draw + * to [0, rng] and nothing is rejected. Fill res directly. */ + err = viRngUniformBits64(VSL_RNG_METHOD_UNIFORMBITS64_STD, + state->stream, len, + (unsigned MKL_INT64 *)res); + assert(err == VSL_STATUS_OK); - while (n_accepted < len) { - npy_intp k = 0; - npy_intp batchSize = len - n_accepted; + DIST_PRAGMA_VECTOR + for (i = 0; i < len; ++i) + res[i] = lo + (res[i] & mask); - err = viRngUniformBits64(VSL_RNG_METHOD_UNIFORM_STD, state->stream, - batchSize, (unsigned MKL_INT64 *)buf); - assert(err == VSL_STATUS_OK); + return; + } - for (k = 0; k < batchSize; ++k) { - npy_uint64 value = buf[k] & mask; - if (value <= rng) { - res[n_accepted++] = lo + value; - } + /* Draw into res and compact in place; n_accepted never runs ahead of i, + * so the store cannot clobber an unread value. Acceptance is > 1/2. */ + err = viRngUniformBits64(VSL_RNG_METHOD_UNIFORMBITS64_STD, + state->stream, len, (unsigned MKL_INT64 *)res); + assert(err == VSL_STATUS_OK); + + for (i = 0; i < len; ++i) { + npy_uint64 value = res[i] & mask; + if (value <= rng) { + res[n_accepted++] = lo + value; } } - mkl_free(buf); + if (n_accepted < len) { + buf = (npy_uint64 *)mkl_malloc( + (len - n_accepted) * sizeof(npy_uint64), 64); + assert(buf != nullptr); + + while (n_accepted < len) { + npy_intp k = 0; + npy_intp batchSize = len - n_accepted; + + err = viRngUniformBits64(VSL_RNG_METHOD_UNIFORMBITS64_STD, + state->stream, batchSize, + (unsigned MKL_INT64 *)buf); + assert(err == VSL_STATUS_OK); + + for (k = 0; k < batchSize; ++k) { + npy_uint64 value = buf[k] & mask; + if (value <= rng) { + res[n_accepted++] = lo + value; + } + } + } + + mkl_free(buf); + } } } diff --git a/mkl_random/tests/test_random.py b/mkl_random/tests/test_random.py index b14fdb3..3025786 100644 --- a/mkl_random/tests/test_random.py +++ b/mkl_random/tests/test_random.py @@ -303,6 +303,37 @@ def test_in_bounds_fuzz(self, randint): assert vals.max() < 2 assert vals.min() >= 0 + @pytest.mark.parametrize( + "high", + [2**31, 2**31 + 1, 2**32, 2**48, 2**63, 2**64 - 1], + ids=["2**31", "2**31+1", "2**32", "2**48", "2**63", "2**64-1"], + ) + @pytest.mark.parametrize("low", [0, 100]) + def test_wide_range_in_bounds(self, low, high): + # Ranges at or above INT_MAX take the masked branch, which + # test_in_bounds_fuzz never reaches (it only uses high <= 16). + for dtype in ("int64", "uint64"): + if low + high > np.iinfo(dtype).max: + continue + vals = rnd.MKLRandomState(1234).randint( + low, high, size=2**20, dtype=dtype + ) + assert vals.min() >= low, f"{dtype}: {vals.min()} < {low}" + assert vals.max() <= high - 1, f"{dtype}: {vals.max()} > {high - 1}" + + @pytest.mark.parametrize( + "dtype", + ["bool", "uint8", "int8", "uint16", "int16", "uint32", "uint64"], + ) + def test_narrow_width_full_range_in_bounds(self, dtype): + # Narrow fills stage in tiles; cover the 4096 tile boundary and beyond. + hi = 2 if dtype == "bool" else int(np.iinfo(dtype).max) + for size in (1, 4095, 4096, 4097, 100000): + vals = rnd.MKLRandomState(7).randint(0, hi, size=size, dtype=dtype) + assert vals.shape == (size,) + assert vals.min() >= 0 + assert vals.max() <= hi - 1 + def test_repeatability(self, randint): import hashlib