From d53e1036e8c6992cd64f5e553a14dcd9c02974ab Mon Sep 17 00:00:00 2001 From: vchamarthi Date: Thu, 3 Sep 2026 22:50:53 -0500 Subject: [PATCH 1/2] fix: derive the 64-bit randint mask from the range maximum --- CHANGELOG.md | 6 + mkl_random/src/mkl_distributions.cpp | 218 +++++++++++---------------- mkl_random/tests/test_random.py | 32 ++++ 3 files changed, 129 insertions(+), 127 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0346ac3b..8af6d913 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,12 +7,18 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 # [dev] (MM/DD/YYYY) ### Added +* Added bound and tile-boundary tests for the integer fills ### 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`; streams for those ranges change +* Sped up `randint` for `bool`, `uint8`, `int8`, `uint16` and `int16`; generated values are unchanged ### 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` +* Fixed the integer fills silently under-filling requests larger than two `MKL_INT_MAX` chunks +* Fixed `multinomial` under-filling large outputs by decrementing its chunk counter by elements instead of draws ## [1.5.0] (08/12/2026) diff --git a/mkl_random/src/mkl_distributions.cpp b/mkl_random/src/mkl_distributions.cpp index 8fa67d96..c6c749f0 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,45 @@ 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) +{ + 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; + int 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 +1785,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 +1794,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 +1808,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 +1817,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 +1831,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 +1840,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 +1854,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 +1863,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 +1877,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 +1892,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 +1946,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 +1983,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 +1999,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 +2010,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 +2037,60 @@ 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 671dcd56..07e1e746 100644 --- a/mkl_random/tests/test_random.py +++ b/mkl_random/tests/test_random.py @@ -275,6 +275,38 @@ def test_randint_in_bounds_fuzz(randint): 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_randint_wide_range_in_bounds(low, high): + # Ranges at or above INT_MAX take the masked branch, which + # test_randint_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_randint_narrow_width_full_range_in_bounds(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_randint_repeatability(randint): import hashlib From 02033af844839d19760c49ecc44968ba2562f39c Mon Sep 17 00:00:00 2001 From: vchamarthi Date: Fri, 4 Sep 2026 09:24:08 -0500 Subject: [PATCH 2/2] fix pr suggestions --- CHANGELOG.md | 12 ++++++------ mkl_random/src/mkl_distributions.cpp | 13 ++++++------- 2 files changed, 12 insertions(+), 13 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8af6d913..0607b857 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,18 +7,18 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 # [dev] (MM/DD/YYYY) ### Added -* Added bound and tile-boundary tests for the integer fills ### 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`; streams for those ranges change -* Sped up `randint` for `bool`, `uint8`, `int8`, `uint16` and `int16`; generated values are unchanged +* 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` -* Fixed the integer fills silently under-filling requests larger than two `MKL_INT_MAX` chunks -* Fixed `multinomial` under-filling large outputs by decrementing its chunk counter by elements instead of draws +* 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 c6c749f0..2b998d8d 100644 --- a/mkl_random/src/mkl_distributions.cpp +++ b/mkl_random/src/mkl_distributions.cpp @@ -1746,14 +1746,15 @@ static void irk_rand_narrow_fill(irk_state *state, 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; - int err = viRngUniform(VSL_RNG_METHOD_UNIFORM_STD, state->stream, - (int)n, tile, lo, hi + 1); + err = viRngUniform(VSL_RNG_METHOD_UNIFORM_STD, state->stream, (int)n, + tile, lo, hi + 1); assert(err == VSL_STATUS_OK); DIST_PRAGMA_VECTOR @@ -2055,8 +2056,7 @@ void irk_rand_uint64_vec(irk_state *state, /* 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); + state->stream, len, (unsigned MKL_INT64 *)res); assert(err == VSL_STATUS_OK); for (i = 0; i < len; ++i) { @@ -2067,9 +2067,8 @@ void irk_rand_uint64_vec(irk_state *state, } if (n_accepted < len) { - buf = (npy_uint64 *)mkl_malloc((len - n_accepted) * - sizeof(npy_uint64), - 64); + buf = (npy_uint64 *)mkl_malloc( + (len - n_accepted) * sizeof(npy_uint64), 64); assert(buf != nullptr); while (n_accepted < len) {