From eb69ce0022af7d51b241909e35a248f1e56df18e Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 16 Oct 2026 05:30:06 -0700 Subject: [PATCH 1/7] Count trials in geometric, as the docstring and NumPy do viRngGeometric returns the number of failures, with support 0, 1, 2, ... and mean 1 / p - 1, and that value was returned unchanged. The docstring documents the NumPy convention instead: "supported on the positive integers, k = 1, 2, ...", a PMF of (1 - p)^(k - 1) * p, and an example answering "How many trials succeeded after a single run?" with p itself. Measured over 200k samples at p = 0.35, the frequency of 1 was p * (1 - p) = 0.2287, and p = 0.3487 was the frequency of 0. This is not an edge case of p == 1, which merely made the mismatch obvious by returning 0 where NumPy returns 1: every p was off by one. Adding the successful trial fixes the whole distribution, so the special case for p == 1 now fills ones. Values that would overflow int are left as they are; they are already meaningless, as MKL writes INT_MIN once 1 / p exceeds the range of its int output. Samples for a given seed therefore grow by one, which the regression test records. The new test pins the convention itself rather than a stream. --- mkl_random/src/mkl_distributions.cpp | 11 +++++++++-- mkl_random/tests/test_random.py | 16 +++++++++++++++- 2 files changed, 24 insertions(+), 3 deletions(-) diff --git a/mkl_random/src/mkl_distributions.cpp b/mkl_random/src/mkl_distributions.cpp index 81c6f0e..7854204 100644 --- a/mkl_random/src/mkl_distributions.cpp +++ b/mkl_random/src/mkl_distributions.cpp @@ -1306,6 +1306,8 @@ void irk_multinomial_vec(irk_state *state, void irk_geometric_vec(irk_state *state, npy_intp len, int *res, const double p) { int err = 0; + npy_intp i = 0, total = len; + int *out = res; if (len < 1) return; @@ -1323,12 +1325,17 @@ void irk_geometric_vec(irk_state *state, npy_intp len, int *res, const double p) err = viRngGeometric(VSL_RNG_METHOD_GEOMETRIC_ICDF, state->stream, len, res, p); assert(err == VSL_STATUS_OK); + + /* MKL counts the failures, the number of trials also counts the + success that ended them; values that would overflow are kept */ + DIST_PRAGMA_VECTOR + for (i = 0; i < total; ++i) + out[i] += (out[i] < INT_MAX) ? 1 : 0; } else { if (p == 1.0) { - npy_intp i; for (i = 0; i < len; ++i) - res[i] = 0; + res[i] = 1; } else { assert(p >= 0.0); diff --git a/mkl_random/tests/test_random.py b/mkl_random/tests/test_random.py index 13eb55b..c2e11ff 100644 --- a/mkl_random/tests/test_random.py +++ b/mkl_random/tests/test_random.py @@ -922,10 +922,24 @@ def test_randomdist_gamma(randomdist): def test_randomdsit_geometric(randomdist): rnd.seed(randomdist.seed, brng=randomdist.brng) actual = rnd.geometric(0.123456789, size=(3, 2)) - desired = np.array([[0, 30], [13, 2], [4, 0]]) + desired = np.array([[1, 31], [14, 3], [5, 1]]) np.testing.assert_array_equal(actual, desired) +def test_geometric_counts_trials(): + # The number of trials, as the docstring and NumPy define it, so the + # support starts at one, P(1) is p and the mean is 1 / p + p = 0.35 + samples = rnd.MKLRandomState(1).geometric(p, 100_000) + assert samples.min() == 1 + np.testing.assert_allclose((samples == 1).mean(), p, atol=0.01) + np.testing.assert_allclose(samples.mean(), 1 / p, rtol=0.01) + + ones = np.ones(5, dtype="int32") + np.testing.assert_array_equal(rnd.geometric(1.0, 5), ones) + np.testing.assert_array_equal(rnd.geometric(np.ones(5)), ones) + + def test_randomdist_gumbel(randomdist): rnd.seed(randomdist.seed, randomdist.brng) actual = rnd.gumbel(loc=0.123456789, scale=2.0, size=(3, 2)) From b3b8f7a76ae00190347b3852dace819197b698db Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 16 Oct 2026 05:32:25 -0700 Subject: [PATCH 2/7] Add gh-196 to changelog --- CHANGELOG.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index b05a04e..78f42a3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -29,8 +29,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 * Sped up `randint` with `array_like` bounds by generating words in cache-sized chunks and by picking the Lemire rejection test by range width [gh-173](https://github.com/IntelPython/mkl_random/pull/173) * Reduced the memory use and sped up integer generation with `WH`, `MCG31`, `R250` and `MRG32K3A`, which lack `viRngUniformBits` support [gh-178](https://github.com/IntelPython/mkl_random/pull/178) * Replaced the `mklrand_py_helper.h` shim with `cpython.bytes` cimports and dropped the no-op `__future__` import, since the project requires Python 3.10+ [gh-192](https://github.com/IntelPython/mkl_random/pull/192) +* `geometric` samples are now one larger than before for every `p`: with a fixed seed they change accordingly. All other streams are unaffected [gh-194](https://github.com/IntelPython/mkl_random/pull/194) ### Fixed +* Fixed `geometric` returning the number of failures before the first success instead of the number of trials, which contradicted its own docstring and NumPy. The support now starts at `1`, `P(1)` is `p` and the mean is `1 / p` [gh-194](https://github.com/IntelPython/mkl_random/pull/194) * Fixed `multinormal_cholesky` and the other functions taking `method` routing integer-like ids such as `np.int64(0)`, `np.int64(2)` or `False` to `BoxMuller`; NumPy integers now select the method they name [gh-189](https://github.com/IntelPython/mkl_random/pull/189) * Fixed `multinormal_cholesky` raising `ZeroDivisionError` for an empty `mean`; it now raises `ValueError` [gh-189](https://github.com/IntelPython/mkl_random/pull/189) * Fixed an out-of-range integer `brng` indexing `brng_list` past its end, which returned uninitialized memory as random values; it now warns and falls back to `MT19937` [gh-177](https://github.com/IntelPython/mkl_random/pull/177) From 643d2a55f035855c24ede3204298738c537e299e Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 16 Oct 2026 06:20:24 -0700 Subject: [PATCH 3/7] Link the changelog entries to this PR gh-194 is the unrelated np.long issue. The entries were written before the PR existed and kept the guessed number, the same slip as in gh-195. --- CHANGELOG.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 78f42a3..1202671 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -29,10 +29,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 * Sped up `randint` with `array_like` bounds by generating words in cache-sized chunks and by picking the Lemire rejection test by range width [gh-173](https://github.com/IntelPython/mkl_random/pull/173) * Reduced the memory use and sped up integer generation with `WH`, `MCG31`, `R250` and `MRG32K3A`, which lack `viRngUniformBits` support [gh-178](https://github.com/IntelPython/mkl_random/pull/178) * Replaced the `mklrand_py_helper.h` shim with `cpython.bytes` cimports and dropped the no-op `__future__` import, since the project requires Python 3.10+ [gh-192](https://github.com/IntelPython/mkl_random/pull/192) -* `geometric` samples are now one larger than before for every `p`: with a fixed seed they change accordingly. All other streams are unaffected [gh-194](https://github.com/IntelPython/mkl_random/pull/194) +* `geometric` samples are now one larger than before for every `p`: with a fixed seed they change accordingly. All other streams are unaffected [gh-196](https://github.com/IntelPython/mkl_random/pull/196) ### Fixed -* Fixed `geometric` returning the number of failures before the first success instead of the number of trials, which contradicted its own docstring and NumPy. The support now starts at `1`, `P(1)` is `p` and the mean is `1 / p` [gh-194](https://github.com/IntelPython/mkl_random/pull/194) +* Fixed `geometric` returning the number of failures before the first success instead of the number of trials, which contradicted its own docstring and NumPy. The support now starts at `1`, `P(1)` is `p` and the mean is `1 / p` [gh-196](https://github.com/IntelPython/mkl_random/pull/196) * Fixed `multinormal_cholesky` and the other functions taking `method` routing integer-like ids such as `np.int64(0)`, `np.int64(2)` or `False` to `BoxMuller`; NumPy integers now select the method they name [gh-189](https://github.com/IntelPython/mkl_random/pull/189) * Fixed `multinormal_cholesky` raising `ZeroDivisionError` for an empty `mean`; it now raises `ValueError` [gh-189](https://github.com/IntelPython/mkl_random/pull/189) * Fixed an out-of-range integer `brng` indexing `brng_list` past its end, which returned uninitialized memory as random values; it now warns and falls back to `MT19937` [gh-177](https://github.com/IntelPython/mkl_random/pull/177) From 2d7ce49d0a4673edad56e77a590a4832d50c5cbd Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 16 Oct 2026 06:42:07 -0700 Subject: [PATCH 4/7] Keep MKL's overflow marker in geometric Per review: the guard never fired. MKL writes INT_MIN on overflow, and INT_MIN < INT_MAX, so the +1 turned it into INT_MIN + 1, which contradicts the comment. Measured over 500k draws at p from 1e-9 down to 1e-12, the only negative value present was INT_MIN + 1 and INT_MAX never appeared; with the sign check the marker now survives unchanged. The comment also explains why total and out exist: the chunk loop advances res and decrements len, so without them the +1 pass would only cover the last chunk for outputs above MKL_INT_MAX, which nothing tests. --- mkl_random/src/mkl_distributions.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/mkl_random/src/mkl_distributions.cpp b/mkl_random/src/mkl_distributions.cpp index 7854204..01adcb2 100644 --- a/mkl_random/src/mkl_distributions.cpp +++ b/mkl_random/src/mkl_distributions.cpp @@ -1306,6 +1306,8 @@ void irk_multinomial_vec(irk_state *state, void irk_geometric_vec(irk_state *state, npy_intp len, int *res, const double p) { int err = 0; + /* res and len advance in the chunk loop, so the +1 pass below + needs the whole range */ npy_intp i = 0, total = len; int *out = res; @@ -1327,10 +1329,11 @@ void irk_geometric_vec(irk_state *state, npy_intp len, int *res, const double p) assert(err == VSL_STATUS_OK); /* MKL counts the failures, the number of trials also counts the - success that ended them; values that would overflow are kept */ + success that ended them; overflowed values, which MKL writes as + negative, are kept */ DIST_PRAGMA_VECTOR for (i = 0; i < total; ++i) - out[i] += (out[i] < INT_MAX) ? 1 : 0; + out[i] += (out[i] >= 0 && out[i] < INT_MAX) ? 1 : 0; } else { if (p == 1.0) { From 4c4556226423bf0a259aac70586c2c84ef4ac6ae Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 16 Oct 2026 06:42:16 -0700 Subject: [PATCH 5/7] Compare dtypes in the geometric ones check Per review: assert_array_equal ignores dtypes by default, so the int32 expectation was not actually checked and an int64 or float result would have passed. strict=True also pins the shape and needs NumPy 1.24, below the 1.26.4 minimum. --- mkl_random/tests/test_random.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mkl_random/tests/test_random.py b/mkl_random/tests/test_random.py index c2e11ff..d487c33 100644 --- a/mkl_random/tests/test_random.py +++ b/mkl_random/tests/test_random.py @@ -936,8 +936,8 @@ def test_geometric_counts_trials(): np.testing.assert_allclose(samples.mean(), 1 / p, rtol=0.01) ones = np.ones(5, dtype="int32") - np.testing.assert_array_equal(rnd.geometric(1.0, 5), ones) - np.testing.assert_array_equal(rnd.geometric(np.ones(5)), ones) + np.testing.assert_array_equal(rnd.geometric(1.0, 5), ones, strict=True) + np.testing.assert_array_equal(rnd.geometric(np.ones(5)), ones, strict=True) def test_randomdist_gumbel(randomdist): From 05fe07fa206ee1f1ee54090999e350baa334cc0b Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 16 Oct 2026 06:42:21 -0700 Subject: [PATCH 6/7] Point the changelog entry at the old values Per review: "with a fixed seed they change accordingly" mostly repeated "one larger", so the entry now says how to get the previous samples back. --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1202671..cdfaa10 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -29,7 +29,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 * Sped up `randint` with `array_like` bounds by generating words in cache-sized chunks and by picking the Lemire rejection test by range width [gh-173](https://github.com/IntelPython/mkl_random/pull/173) * Reduced the memory use and sped up integer generation with `WH`, `MCG31`, `R250` and `MRG32K3A`, which lack `viRngUniformBits` support [gh-178](https://github.com/IntelPython/mkl_random/pull/178) * Replaced the `mklrand_py_helper.h` shim with `cpython.bytes` cimports and dropped the no-op `__future__` import, since the project requires Python 3.10+ [gh-192](https://github.com/IntelPython/mkl_random/pull/192) -* `geometric` samples are now one larger than before for every `p`: with a fixed seed they change accordingly. All other streams are unaffected [gh-196](https://github.com/IntelPython/mkl_random/pull/196) +* `geometric` samples are now one larger than before for every `p`, so results for a fixed seed change; `geometric(p) - 1` reproduces the old values. All other streams are unaffected [gh-196](https://github.com/IntelPython/mkl_random/pull/196) ### Fixed * Fixed `geometric` returning the number of failures before the first success instead of the number of trials, which contradicted its own docstring and NumPy. The support now starts at `1`, `P(1)` is `p` and the mean is `1 / p` [gh-196](https://github.com/IntelPython/mkl_random/pull/196) From 485838294f8a15b026b3debf37f1b9e8ac6550e8 Mon Sep 17 00:00:00 2001 From: Vladislav Perevezentsev Date: Fri, 16 Oct 2026 09:32:38 -0700 Subject: [PATCH 7/7] Pin MKL's overflow marker in the geometric test Per review: the +1 pass skips negatives so the marker survives, but nothing checked it. At p = 1e-300 every draw overflows, so the assertion is exact rather than statistical; measured as all INT_MIN for p of 1e-100, 1e-300 and the smallest subnormal, on both the scalar and the array path. --- mkl_random/tests/test_random.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/mkl_random/tests/test_random.py b/mkl_random/tests/test_random.py index d487c33..3e48f7a 100644 --- a/mkl_random/tests/test_random.py +++ b/mkl_random/tests/test_random.py @@ -939,6 +939,9 @@ def test_geometric_counts_trials(): np.testing.assert_array_equal(rnd.geometric(1.0, 5), ones, strict=True) np.testing.assert_array_equal(rnd.geometric(np.ones(5)), ones, strict=True) + # MKL marks draws that overflow int with INT_MIN, which must survive the +1 + assert (rnd.geometric(1e-300, 3) == np.iinfo(np.int32).min).all() + def test_randomdist_gumbel(randomdist): rnd.seed(randomdist.seed, randomdist.brng)