From f9244770b1991406e585094f8f50007efb411e04 Mon Sep 17 00:00:00 2001 From: Python Infra CI Date: Mon, 5 Oct 2026 18:03:42 +0000 Subject: [PATCH 1/3] fix: chunk irk_multinormal_vec_ICDF draws at MKL_INT_MAX like its BM1/BM --- mkl_random/src/mkl_distributions.cpp | 13 ++++++++++ mkl_random/tests/test_random.py | 38 ++++++++++++++++++++++++++++ 2 files changed, 51 insertions(+) diff --git a/mkl_random/src/mkl_distributions.cpp b/mkl_random/src/mkl_distributions.cpp index 1427610..29a0ff7 100644 --- a/mkl_random/src/mkl_distributions.cpp +++ b/mkl_random/src/mkl_distributions.cpp @@ -2287,6 +2287,19 @@ void irk_multinormal_vec_ICDF(irk_state *state, int err = 0; const MKL_INT storage_mode = cholesky_storage_flags[storage_flag]; + if (len < 1) + return; + + while (len > MKL_INT_MAX) { + err = + vdRngGaussianMV(VSL_RNG_METHOD_GAUSSIANMV_ICDF, state->stream, + MKL_INT_MAX, res, dim, storage_mode, mean_vec, ch); + assert(err == VSL_STATUS_OK); + + res += MKL_INT_MAX * dim; + len -= MKL_INT_MAX; + } + err = vdRngGaussianMV(VSL_RNG_METHOD_GAUSSIANMV_ICDF, state->stream, len, res, dim, storage_mode, mean_vec, ch); assert(err == VSL_STATUS_OK); diff --git a/mkl_random/tests/test_random.py b/mkl_random/tests/test_random.py index 89050da..51af2c0 100644 --- a/mkl_random/tests/test_random.py +++ b/mkl_random/tests/test_random.py @@ -1119,6 +1119,44 @@ def test_randomdist_multinormal_cholesky(randomdist): np.testing.assert_allclose(actual, desired, atol=1e-10, rtol=1e-10) +def test_multinormal_cholesky_icdf_matches_split_calls(): + # irk_multinormal_vec_ICDF had no MKL_INT_MAX chunking loop, unlike its + # BM1/BM2 siblings (see mkl_distributions.cpp). The overflow itself only + # shows up for draw counts above MKL_INT_MAX (~2**31), far too large to + # allocate in a test, so this instead checks the invariant the chunking + # loop relies on: splitting one request into two consecutive calls that + # share the stream must equal one call for the combined count. The fix + # does not change this result; it would only matter above MKL_INT_MAX. + mean = np.array([0.1, -0.2]) + chol_mat = np.array([[1.0, 0.0], [-0.5, 1.0]]) + n1, n2 = 17, 29 + + whole = rnd.MKLRandomState(123).multinormal_cholesky( + mean, chol_mat, size=n1 + n2, method="ICDF" + ) + + split_state = rnd.MKLRandomState(123) + first = split_state.multinormal_cholesky( + mean, chol_mat, size=n1, method="ICDF" + ) + second = split_state.multinormal_cholesky( + mean, chol_mat, size=n2, method="ICDF" + ) + + np.testing.assert_array_equal(whole[:n1], first) + np.testing.assert_array_equal(whole[n1:], second) + + +def test_multinormal_cholesky_icdf_empty_size(): + # len == 0 must be a no-op rather than issuing a 0-count MKL call. + mean = np.array([0.1, -0.2]) + chol_mat = np.array([[1.0, 0.0], [-0.5, 1.0]]) + out = rnd.MKLRandomState(123).multinormal_cholesky( + mean, chol_mat, size=0, method="ICDF" + ) + assert out.shape == (0, 2) + + def test_randomdist_negative_binomial(randomdist): rnd.seed(randomdist.seed, brng=randomdist.brng) actual = rnd.negative_binomial(n=100, p=0.12345, size=(3, 2)) From 4df72afe835490d3c203f4da83c2c00ac7142af2 Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Tue, 6 Oct 2026 13:38:52 +0200 Subject: [PATCH 2/3] test: tidy multinormal_cholesky ICDF tests and add changelog entry - Shorten the split-call test comment and compare with assert_allclose, since bit-exact continuation across calls is an MKL implementation detail that does not hold for every BRNG. - Reword the empty-size test: a 0-count MKL call was already a no-op, so check instead that an empty draw returns an empty array and leaves the stream untouched, for all three methods and for size=0 and (3, 0). - Add a Fixed entry to CHANGELOG.md. --- CHANGELOG.md | 1 + mkl_random/tests/test_random.py | 28 ++++++++++++++-------------- 2 files changed, 15 insertions(+), 14 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ecf2d23..2e53a54 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -39,6 +39,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 * Fixed `randint` with scalar bounds above `INT_MAX`, `tomaxint` and `bytes` returning uninitialized data for `WH`, `MCG31`, `R250` and `MRG32K3A`, which lack `viRngUniformBits` support [gh-175](https://github.com/IntelPython/mkl_random/pull/175) * Fixed 64-bit integer generation returning zeros or crashing with `PHILOX4X32X10` and `ARS5` for requests of `2**30` elements or more [gh-175](https://github.com/IntelPython/mkl_random/pull/175) * Fixed `randint_untyped` raising `OverflowError` for bounds outside the C `long` range, e.g. `2**40` on Windows [gh-176](https://github.com/IntelPython/mkl_random/pull/176) +* Fixed `multinormal_cholesky` with `method="ICDF"` silently leaving the output uninitialized or partially filled for requests of more than `MKL_INT_MAX` samples [gh-187](https://github.com/IntelPython/mkl_random/pull/187) ### Removed * Removed the `python-gil` constraint from the conda recipes, which pinned `mkl_random` to GIL-enabled Python 3.14 builds [gh-159](https://github.com/IntelPython/mkl_random/pull/159) diff --git a/mkl_random/tests/test_random.py b/mkl_random/tests/test_random.py index 51af2c0..e94c029 100644 --- a/mkl_random/tests/test_random.py +++ b/mkl_random/tests/test_random.py @@ -1120,13 +1120,8 @@ def test_randomdist_multinormal_cholesky(randomdist): def test_multinormal_cholesky_icdf_matches_split_calls(): - # irk_multinormal_vec_ICDF had no MKL_INT_MAX chunking loop, unlike its - # BM1/BM2 siblings (see mkl_distributions.cpp). The overflow itself only - # shows up for draw counts above MKL_INT_MAX (~2**31), far too large to - # allocate in a test, so this instead checks the invariant the chunking - # loop relies on: splitting one request into two consecutive calls that - # share the stream must equal one call for the combined count. The fix - # does not change this result; it would only matter above MKL_INT_MAX. + # Results must not depend on how a draw is split across calls; the C + # layer relies on this when it chunks requests at MKL_INT_MAX. mean = np.array([0.1, -0.2]) chol_mat = np.array([[1.0, 0.0], [-0.5, 1.0]]) n1, n2 = 17, 29 @@ -1143,18 +1138,23 @@ def test_multinormal_cholesky_icdf_matches_split_calls(): mean, chol_mat, size=n2, method="ICDF" ) - np.testing.assert_array_equal(whole[:n1], first) - np.testing.assert_array_equal(whole[n1:], second) + np.testing.assert_allclose(whole[:n1], first, rtol=1e-12, atol=1e-12) + np.testing.assert_allclose(whole[n1:], second, rtol=1e-12, atol=1e-12) -def test_multinormal_cholesky_icdf_empty_size(): - # len == 0 must be a no-op rather than issuing a 0-count MKL call. +@pytest.mark.parametrize("method", ["ICDF", "BoxMuller", "BoxMuller2"]) +@pytest.mark.parametrize("size,shape", [(0, (0, 2)), ((3, 0), (3, 0, 2))]) +def test_multinormal_cholesky_empty_size(method, size, shape): + # An empty draw returns an empty array and leaves the stream untouched. mean = np.array([0.1, -0.2]) chol_mat = np.array([[1.0, 0.0], [-0.5, 1.0]]) - out = rnd.MKLRandomState(123).multinormal_cholesky( - mean, chol_mat, size=0, method="ICDF" + reference = rnd.MKLRandomState(123) + state = rnd.MKLRandomState(123) + out = state.multinormal_cholesky(mean, chol_mat, size=size, method=method) + assert out.shape == shape + np.testing.assert_array_equal( + state.random_sample(32), reference.random_sample(32) ) - assert out.shape == (0, 2) def test_randomdist_negative_binomial(randomdist): From 42a8f947f1a2a61841a8ebefaab5f7aab85be613 Mon Sep 17 00:00:00 2001 From: Anton Volkov Date: Tue, 6 Oct 2026 19:54:16 +0200 Subject: [PATCH 3/3] fix: chunk multinormal and multinomial draws so each oneMKL call fits in MKL_INT oneMKL's 32-bit interface mishandles a single vdRngGaussianMV call once n * dim exceeds MKL_INT_MAX: it returns VSL_STATUS_OK but generates nothing (or only a wrapped-around prefix). Chunking at MKL_INT_MAX vectors, as irk_multinormal_vec_BM1/BM2 did, is therefore not enough for dim >= 2. - Chunk irk_multinormal_vec_ICDF/BM1/BM2 at MKL_INT_MAX / dim vectors and advance the output pointer by max_len * dim doubles. - Chunk irk_multinomial_vec at MKL_INT_MAX / k as a defensive measure. - Run the split-call test for all three methods. - Reword the CHANGELOG entry to cover the general n * dim condition. --- CHANGELOG.md | 2 +- mkl_random/src/mkl_distributions.cpp | 51 ++++++++++++++++++---------- mkl_random/tests/test_random.py | 11 +++--- 3 files changed, 40 insertions(+), 24 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2e53a54..a7c982f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -39,7 +39,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 * Fixed `randint` with scalar bounds above `INT_MAX`, `tomaxint` and `bytes` returning uninitialized data for `WH`, `MCG31`, `R250` and `MRG32K3A`, which lack `viRngUniformBits` support [gh-175](https://github.com/IntelPython/mkl_random/pull/175) * Fixed 64-bit integer generation returning zeros or crashing with `PHILOX4X32X10` and `ARS5` for requests of `2**30` elements or more [gh-175](https://github.com/IntelPython/mkl_random/pull/175) * Fixed `randint_untyped` raising `OverflowError` for bounds outside the C `long` range, e.g. `2**40` on Windows [gh-176](https://github.com/IntelPython/mkl_random/pull/176) -* Fixed `multinormal_cholesky` with `method="ICDF"` silently leaving the output uninitialized or partially filled for requests of more than `MKL_INT_MAX` samples [gh-187](https://github.com/IntelPython/mkl_random/pull/187) +* Fixed `multinormal_cholesky` silently returning wrong output, e.g. every sample equal to `mean`, when the number of samples times the dimension exceeds `MKL_INT_MAX` [gh-187](https://github.com/IntelPython/mkl_random/pull/187) ### Removed * Removed the `python-gil` constraint from the conda recipes, which pinned `mkl_random` to GIL-enabled Python 3.14 builds [gh-159](https://github.com/IntelPython/mkl_random/pull/159) diff --git a/mkl_random/src/mkl_distributions.cpp b/mkl_random/src/mkl_distributions.cpp index 29a0ff7..9320ad7 100644 --- a/mkl_random/src/mkl_distributions.cpp +++ b/mkl_random/src/mkl_distributions.cpp @@ -1284,13 +1284,17 @@ void irk_multinomial_vec(irk_state *state, memset(res, 0, len * k * sizeof(int)); } else { - while (len > MKL_INT_MAX) { + /* Keep len * k within MKL_INT, since oneMKL indexes the output + in MKL_INT. */ + const npy_intp max_len = MKL_INT_MAX / k; + + while (len > max_len) { err = viRngMultinomial(VSL_RNG_METHOD_MULTINOMIAL_MULTPOISSON, - state->stream, MKL_INT_MAX, res, n, k, pvec); + state->stream, max_len, res, n, k, pvec); assert(err == VSL_STATUS_OK); /* len counts draws, res counts ints. */ - res += k * MKL_INT_MAX; - len -= MKL_INT_MAX; + res += k * max_len; + len -= max_len; } err = viRngMultinomial(VSL_RNG_METHOD_MULTINOMIAL_MULTPOISSON, @@ -2290,14 +2294,17 @@ void irk_multinormal_vec_ICDF(irk_state *state, if (len < 1) return; - while (len > MKL_INT_MAX) { - err = - vdRngGaussianMV(VSL_RNG_METHOD_GAUSSIANMV_ICDF, state->stream, - MKL_INT_MAX, res, dim, storage_mode, mean_vec, ch); + /* oneMKL fills len * dim doubles per call, and that count must fit in + MKL_INT. */ + const npy_intp max_len = MKL_INT_MAX / dim; + + while (len > max_len) { + err = vdRngGaussianMV(VSL_RNG_METHOD_GAUSSIANMV_ICDF, state->stream, + max_len, res, dim, storage_mode, mean_vec, ch); assert(err == VSL_STATUS_OK); - res += MKL_INT_MAX * dim; - len -= MKL_INT_MAX; + res += max_len * dim; + len -= max_len; } err = vdRngGaussianMV(VSL_RNG_METHOD_GAUSSIANMV_ICDF, state->stream, len, @@ -2319,14 +2326,18 @@ void irk_multinormal_vec_BM1(irk_state *state, if (len < 1) return; - while (len > MKL_INT_MAX) { + /* oneMKL fills len * dim doubles per call, and that count must fit in + MKL_INT. */ + const npy_intp max_len = MKL_INT_MAX / dim; + + while (len > max_len) { err = vdRngGaussianMV(VSL_RNG_METHOD_GAUSSIANMV_BOXMULLER, state->stream, - MKL_INT_MAX, res, dim, storage_mode, mean_vec, ch); + max_len, res, dim, storage_mode, mean_vec, ch); assert(err == VSL_STATUS_OK); - res += MKL_INT_MAX * dim; - len -= MKL_INT_MAX; + res += max_len * dim; + len -= max_len; } err = vdRngGaussianMV(VSL_RNG_METHOD_GAUSSIANMV_BOXMULLER, state->stream, @@ -2348,14 +2359,18 @@ void irk_multinormal_vec_BM2(irk_state *state, if (len < 1) return; - while (len > MKL_INT_MAX) { + /* oneMKL fills len * dim doubles per call, and that count must fit in + MKL_INT. */ + const npy_intp max_len = MKL_INT_MAX / dim; + + while (len > max_len) { err = vdRngGaussianMV(VSL_RNG_METHOD_GAUSSIANMV_BOXMULLER2, state->stream, - MKL_INT_MAX, res, dim, storage_mode, mean_vec, ch); + max_len, res, dim, storage_mode, mean_vec, ch); assert(err == VSL_STATUS_OK); - res += MKL_INT_MAX * dim; - len -= MKL_INT_MAX; + res += max_len * dim; + len -= max_len; } err = vdRngGaussianMV(VSL_RNG_METHOD_GAUSSIANMV_BOXMULLER2, state->stream, diff --git a/mkl_random/tests/test_random.py b/mkl_random/tests/test_random.py index e94c029..d508f4e 100644 --- a/mkl_random/tests/test_random.py +++ b/mkl_random/tests/test_random.py @@ -1119,23 +1119,24 @@ def test_randomdist_multinormal_cholesky(randomdist): np.testing.assert_allclose(actual, desired, atol=1e-10, rtol=1e-10) -def test_multinormal_cholesky_icdf_matches_split_calls(): +@pytest.mark.parametrize("method", ["ICDF", "BoxMuller", "BoxMuller2"]) +def test_multinormal_cholesky_matches_split_calls(method): # Results must not depend on how a draw is split across calls; the C - # layer relies on this when it chunks requests at MKL_INT_MAX. + # layer relies on this when it splits large requests into chunks. mean = np.array([0.1, -0.2]) chol_mat = np.array([[1.0, 0.0], [-0.5, 1.0]]) n1, n2 = 17, 29 whole = rnd.MKLRandomState(123).multinormal_cholesky( - mean, chol_mat, size=n1 + n2, method="ICDF" + mean, chol_mat, size=n1 + n2, method=method ) split_state = rnd.MKLRandomState(123) first = split_state.multinormal_cholesky( - mean, chol_mat, size=n1, method="ICDF" + mean, chol_mat, size=n1, method=method ) second = split_state.multinormal_cholesky( - mean, chol_mat, size=n2, method="ICDF" + mean, chol_mat, size=n2, method=method ) np.testing.assert_allclose(whole[:n1], first, rtol=1e-12, atol=1e-12)