Skip to content

fix: 64-bit randint could return high; correct the integer fills #172

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
vchamarthi wants to merge 4 commits into
base: master
Choose a base branch
from
+127 −127
Open

fix: 64-bit randint could return high; correct the integer fills #172

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
d53e103
fix: derive the 64-bit randint mask from the range maximum
vchamarthi Sep 4, 2026
02033af
fix pr suggestions
vchamarthi Sep 4, 2026
adc5356
Merge remote-tracking branch 'upstream/master' into fix/randint-integ…
vchamarthi Sep 4, 2026
519be4d
Merge branch 'master' into fix/randint-integer-fills
ndgrigorian Sep 9, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 6 additions & 0 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down
217 changes: 90 additions & 127 deletions mkl_random/src/mkl_distributions.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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 <typename T>
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;
Comment thread
vchamarthi marked this conversation as resolved.
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)
Expand All @@ -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<npy_bool>(state, len, res, (int)lo, (int)hi);
}

void irk_rand_uint8_vec(irk_state *state,
Expand All @@ -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)
Expand All @@ -1808,18 +1809,7 @@ void irk_rand_uint8_vec(irk_state *state,
}

assert(lo < hi);
Comment on lines 1803 to 1811

@vlad-perevezentsev vlad-perevezentsev Sep 14, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

looks identical in all five narrow fills so they can move into irk_rand_narrow_fill
note:
lo is an int there so the fill needs a cast

 if (lo == hi) {
    DIST_PRAGMA_VECTOR
    for (i = 0; i < len; ++i)
        res[i] = (T)lo;

    return;
}

assert(lo < hi);

@vlad-perevezentsev vlad-perevezentsev Sep 14, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Then if (len < 1) above on 1800 line looks redundant because the constant fill and while (len > 0)already do nothing forlen <= 0`

@vlad-perevezentsev vlad-perevezentsev Sep 14, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For irk_rand_bool_vec the remaining assert((lo == 0) && (hi == 1)) can either stay or go.
Once lo == hi is handled inside the helper assert(lo < hi) already implies it for npy_bool
If you keep it it has to become assert((lo == hi) || ((lo == 0) && (hi == 1))) because it is now reachable for lo == hi as well

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<npy_uint8>(state, len, res, (int)lo, (int)hi);
}

void irk_rand_int8_vec(irk_state *state,
Expand All @@ -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)
Expand All @@ -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<npy_int8>(state, len, res, (int)lo, (int)hi);
}

void irk_rand_uint16_vec(irk_state *state,
Expand All @@ -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)
Expand All @@ -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<npy_uint16>(state, len, res, (int)lo, (int)hi);
}

void irk_rand_int16_vec(irk_state *state,
Expand All @@ -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)
Expand All @@ -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<npy_int16>(state, len, res, (int)lo, (int)hi);
}

void irk_rand_uint32_vec(irk_state *state,
Expand All @@ -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;
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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);
}
}
}

Expand Down
Loading
Loading