From 1bd97cf3e75911422bd376119d3cf1258fc56428 Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Fri, 16 Oct 2026 05:17:29 -0700
Subject: [PATCH 01/10] Return zeros for poisson with lam == 0 and reject NaN
 lam

viRngPoisson requires lam > 0, so lam == 0 made it fail with "Parameter 5
was incorrect on entry to viRngPoisson" and leave the output buffer
untouched, which was then handed back as random values. The assert on the
MKL status code cannot catch that: meson.build defines NDEBUG, so every
assert in mkl_distributions.cpp is a no-op in shipped builds. The scalar
path hit this directly, and so did array-valued lam with method="PTPE",
which calls the scalar kernel per element; only the default POISNORM
array path was correct, since viRngPoissonV accepts a zero lambda.

The distribution is degenerate at lam == 0, so the kernels now fill zeros
without touching the stream, which is what the array path and NumPy
already returned. Streams for a positive lam are unaffected.

NaN lam passed validation as well, because `lam < 0` is false for NaN,
and MKL then returned arbitrary values. It now raises ValueError with
NumPy's wording. The isnan check runs first, as the comparison ufuncs
warn on NaN input.
---
 mkl_random/mklrand.pyx               | 11 ++++++-----
 mkl_random/src/mkl_distributions.cpp | 16 ++++++++++++++++
 mkl_random/tests/test_random.py      | 15 +++++++++++++++
 3 files changed, 37 insertions(+), 5 deletions(-)

diff --git a/mkl_random/mklrand.pyx b/mkl_random/mklrand.pyx
index 7f55da6..bc8ea88 100644
--- a/mkl_random/mklrand.pyx
+++ b/mkl_random/mklrand.pyx
@@ -43,6 +43,7 @@ from cpython.bytes cimport (
     PyBytes_Check,
     PyBytes_FromStringAndSize,
 )
+from libc.math cimport isnan
 from libc.stdint cimport int64_t
 from libc.string cimport memcpy, memset
 
@@ -6052,9 +6053,9 @@ cdef class _MKLRandomState:
         if cnp.PyArray_NDIM(olam) == 0:
             flam = PyFloat_AsDouble(lam)
 
-            if lam < 0:
-                raise ValueError("lam < 0")
-            if lam > self._poisson_lam_max:
+            if isnan(flam) or flam < 0:
+                raise ValueError("lam < 0 or lam is NaN")
+            if flam > self._poisson_lam_max:
                 raise ValueError("lam value too large")
             method = choose_method(
                 method, [POISNORM, PTPE], _method_alias_dict_poisson
@@ -6076,8 +6077,8 @@ cdef class _MKLRandomState:
                     self.lock
                 )
 
-        if np.any(np.less(olam, 0)):
-            raise ValueError("lam < 0")
+        if np.any(np.isnan(olam)) or np.any(np.less(olam, 0)):
+            raise ValueError("lam < 0 or lam is NaN")
         if np.any(np.greater(olam, self._poisson_lam_max)):
             raise ValueError("lam value too large.")
         method = choose_method(
diff --git a/mkl_random/src/mkl_distributions.cpp b/mkl_random/src/mkl_distributions.cpp
index 81c6f0e..c8ad136 100644
--- a/mkl_random/src/mkl_distributions.cpp
+++ b/mkl_random/src/mkl_distributions.cpp
@@ -1399,6 +1399,14 @@ void irk_poisson_vec_PTPE(irk_state *state,
     if (len < 1)
         return;
 
+    /* MKL requires lambda > 0, while lambda == 0 always gives zero events */
+    if (lambda == 0.0) {
+        npy_intp i;
+        for (i = 0; i < len; ++i)
+            res[i] = 0;
+        return;
+    }
+
     while (len > MKL_INT_MAX) {
         err = viRngPoisson(VSL_RNG_METHOD_POISSON_PTPE, state->stream,
                            MKL_INT_MAX, res, lambda);
@@ -1423,6 +1431,14 @@ void irk_poisson_vec_POISNORM(irk_state *state,
     if (len < 1)
         return;
 
+    /* MKL requires lambda > 0, while lambda == 0 always gives zero events */
+    if (lambda == 0.0) {
+        npy_intp i;
+        for (i = 0; i < len; ++i)
+            res[i] = 0;
+        return;
+    }
+
     while (len > MKL_INT_MAX) {
         err = viRngPoisson(VSL_RNG_METHOD_POISSON_POISNORM, state->stream,
                            MKL_INT_MAX, res, lambda);
diff --git a/mkl_random/tests/test_random.py b/mkl_random/tests/test_random.py
index 13eb55b..4a1d59d 100644
--- a/mkl_random/tests/test_random.py
+++ b/mkl_random/tests/test_random.py
@@ -1282,6 +1282,21 @@ def test_poisson_exceptions():
     assert_raises(ValueError, rnd.poisson, [lamneg] * 10)
     assert_raises(ValueError, rnd.poisson, lambig)
     assert_raises(ValueError, rnd.poisson, [lambig] * 10)
+    assert_raises(ValueError, rnd.poisson, np.nan)
+    assert_raises(ValueError, rnd.poisson, [np.nan] * 10)
+
+
+@pytest.mark.parametrize("method", ["POISNORM", "PTPE"])
+def test_poisson_zero_lam(method):
+    # MKL requires lam > 0, while lam == 0 always gives zero events
+    zeros = np.zeros(5, dtype="int32")
+    assert rnd.poisson(0.0, method=method) == 0
+    np.testing.assert_array_equal(rnd.poisson(0.0, 5, method=method), zeros)
+    np.testing.assert_array_equal(
+        rnd.poisson(np.zeros(5), method=method), zeros
+    )
+    mixed = rnd.poisson(np.array([0.0, 5.0, 0.0]), method=method)
+    np.testing.assert_array_equal(mixed[::2], zeros[:2])
 
 
 def test_randomdist_power(randomdist):

From 54b53ae4c7d006a9631a7dc6f926ed17eb9537b7 Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Fri, 16 Oct 2026 05:18:04 -0700
Subject: [PATCH 02/10] Return zeros for negative_binomial with p == 1 and
 reject p == 0

viRngNegbinomial requires 0 < p < 1, and both endpoints reached it. p == 1
only reported "INTERNAL ERROR: Condition -6", and p == 0 passed a check
that tested `p < 0`, so in both cases the untouched output buffer was
handed back as random values, on the scalar and on the array path alike.
The assert on the MKL status code is a no-op, as meson.build defines
NDEBUG.

p == 1 means no trial ever fails, so the number of failures is zero and
the kernel now fills zeros without touching the stream, as NumPy returns.
p == 0 has no finite answer, NumPy's own Generator rejects it, and it now
raises ValueError, which makes the documented range of p (0, 1]. NaN was
unchecked for both n and p, so it is rejected together with them.
---
 mkl_random/mklrand.pyx               | 28 ++++++++++++++--------------
 mkl_random/src/mkl_distributions.cpp |  8 ++++++++
 mkl_random/tests/test_random.py      | 19 +++++++++++++++++++
 3 files changed, 41 insertions(+), 14 deletions(-)

diff --git a/mkl_random/mklrand.pyx b/mkl_random/mklrand.pyx
index bc8ea88..63bf65b 100644
--- a/mkl_random/mklrand.pyx
+++ b/mkl_random/mklrand.pyx
@@ -5872,14 +5872,14 @@ cdef class _MKLRandomState:
 
         Samples are drawn from a negative binomial distribution with specified
         parameters, `n` trials and `p` probability of success where `n` is an
-        integer > 0 and `p` is in the interval [0, 1].
+        integer > 0 and `p` is in the interval (0, 1].
 
         Parameters
         ----------
         n : int
             Parameter, > 0.
         p : float
-            Parameter, >= 0 and <=1.
+            Parameter, > 0 and <= 1.
         size : int or tuple of ints, optional
             Output shape.  If the given shape is, e.g., ``(m, n, k)``, then
             ``m * n * k`` samples are drawn.  Default is None, in which case a
@@ -5950,12 +5950,10 @@ cdef class _MKLRandomState:
             fp = PyFloat_AsDouble(p)
             fn = PyFloat_AsDouble(n)
 
-            if fn <= 0:
-                raise ValueError("n <= 0")
-            if fp < 0:
-                raise ValueError("p < 0")
-            elif fp > 1:
-                raise ValueError("p > 1")
+            if isnan(fn) or fn <= 0:
+                raise ValueError("n <= 0 or n is NaN")
+            if isnan(fp) or fp <= 0 or fp > 1:
+                raise ValueError("p <= 0, p > 1 or p is NaN")
             return vec_discdd_array_sc(
                 self.internal_state,
                 irk_negbinomial_vec,
@@ -5965,12 +5963,14 @@ cdef class _MKLRandomState:
                 self.lock
             )
 
-        if np.any(np.less_equal(n, 0)):
-            raise ValueError("n <= 0")
-        if np.any(np.less(p, 0)):
-            raise ValueError("p < 0")
-        if np.any(np.greater(p, 1)):
-            raise ValueError("p > 1")
+        if np.any(np.isnan(on)) or np.any(np.less_equal(on, 0)):
+            raise ValueError("n <= 0 or n is NaN")
+        if (
+            np.any(np.isnan(op))
+            or np.any(np.less_equal(op, 0))
+            or np.any(np.greater(op, 1))
+        ):
+            raise ValueError("p <= 0, p > 1 or p is NaN")
         return vec_discdd_array(
             self.internal_state, irk_negbinomial_vec, size, on, op, self.lock
         )
diff --git a/mkl_random/src/mkl_distributions.cpp b/mkl_random/src/mkl_distributions.cpp
index c8ad136..52fcdc7 100644
--- a/mkl_random/src/mkl_distributions.cpp
+++ b/mkl_random/src/mkl_distributions.cpp
@@ -1348,6 +1348,14 @@ void irk_negbinomial_vec(irk_state *state,
     if (len < 1)
         return;
 
+    /* MKL requires p < 1, while p == 1 means every trial succeeds */
+    if (p == 1.0) {
+        npy_intp i;
+        for (i = 0; i < len; ++i)
+            res[i] = 0;
+        return;
+    }
+
     while (len > MKL_INT_MAX) {
         err = viRngNegbinomial(VSL_RNG_METHOD_NEGBINOMIAL_NBAR, state->stream,
                                MKL_INT_MAX, res, a, p);
diff --git a/mkl_random/tests/test_random.py b/mkl_random/tests/test_random.py
index 4a1d59d..7f3b2b0 100644
--- a/mkl_random/tests/test_random.py
+++ b/mkl_random/tests/test_random.py
@@ -1169,6 +1169,25 @@ def test_randomdist_negative_binomial(randomdist):
     np.testing.assert_array_equal(actual, desired)
 
 
+def test_negative_binomial_p_is_one():
+    # MKL requires p < 1, while p == 1 means no trial ever fails
+    zeros = np.zeros(5, dtype="int32")
+    assert rnd.negative_binomial(5, 1.0) == 0
+    np.testing.assert_array_equal(rnd.negative_binomial(5, 1.0, 5), zeros)
+    np.testing.assert_array_equal(
+        rnd.negative_binomial(np.full(5, 5.0), 1.0), zeros
+    )
+
+
+def test_negative_binomial_exceptions():
+    for p in (0.0, -0.5, 1.5, np.nan):
+        assert_raises(ValueError, rnd.negative_binomial, 5, p)
+        assert_raises(ValueError, rnd.negative_binomial, 5, [p] * 5)
+    for n in (0, -1, np.nan):
+        assert_raises(ValueError, rnd.negative_binomial, n, 0.5)
+        assert_raises(ValueError, rnd.negative_binomial, [n] * 5, 0.5)
+
+
 def test_randomdist_noncentral_chisquare(randomdist):
     rnd.seed(randomdist.seed, brng=randomdist.brng)
     actual = rnd.noncentral_chisquare(df=5, nonc=5, size=(3, 2))

From d73be2f3ff20c43a8ba85e6890af56ca52d68565 Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Fri, 16 Oct 2026 05:18:32 -0700
Subject: [PATCH 03/10] Reject NaN p in geometric and logseries

Neither validated NaN: `p <= 0.0`, `p > 1.0` and `p >= 1.0` are all false
for NaN, so it reached the kernels. geometric then wrote nothing and
handed back uninitialized memory as random values, since viRngGeometric
is only called for 0 < p < 1; logseries silently returned ones. NumPy
raises ValueError for both.

The checks are merged into one condition per path, with isnan first, as
the comparison ufuncs warn on NaN input.
---
 mkl_random/mklrand.pyx          | 32 ++++++++++++++++----------------
 mkl_random/tests/test_random.py | 12 ++++++++++++
 2 files changed, 28 insertions(+), 16 deletions(-)

diff --git a/mkl_random/mklrand.pyx b/mkl_random/mklrand.pyx
index 63bf65b..49e0106 100644
--- a/mkl_random/mklrand.pyx
+++ b/mkl_random/mklrand.pyx
@@ -6252,18 +6252,18 @@ cdef class _MKLRandomState:
         if cnp.PyArray_NDIM(op) == 0:
             fp = PyFloat_AsDouble(p)
 
-            if fp <= 0.0:
-                raise ValueError("p <= 0.0")
-            if fp > 1.0:
-                raise ValueError("p > 1.0")
+            if isnan(fp) or fp <= 0.0 or fp > 1.0:
+                raise ValueError("p <= 0, p > 1 or p is NaN")
             return vec_discd_array_sc(
                 self.internal_state, irk_geometric_vec, size, fp, self.lock
             )
 
-        if np.any(np.less_equal(op, 0.0)):
-            raise ValueError("p <= 0.0")
-        if np.any(np.greater(op, 1.0)):
-            raise ValueError("p > 1.0")
+        if (
+            np.any(np.isnan(op))
+            or np.any(np.less_equal(op, 0.0))
+            or np.any(np.greater(op, 1.0))
+        ):
+            raise ValueError("p <= 0, p > 1 or p is NaN")
         return vec_discd_array(
             self.internal_state, irk_geometric_vec, size, op, self.lock
         )
@@ -6514,18 +6514,18 @@ cdef class _MKLRandomState:
         if cnp.PyArray_NDIM(op) == 0:
             fp = PyFloat_AsDouble(p)
 
-            if fp <= 0.0:
-                raise ValueError("p <= 0.0")
-            if fp >= 1.0:
-                raise ValueError("p >= 1.0")
+            if isnan(fp) or fp <= 0.0 or fp >= 1.0:
+                raise ValueError("p <= 0, p >= 1 or p is NaN")
             return vec_discd_array_sc(
                 self.internal_state, irk_logseries_vec, size, fp, self.lock
             )
 
-        if np.any(np.less_equal(op, 0.0)):
-            raise ValueError("p <= 0.0")
-        if np.any(np.greater_equal(op, 1.0)):
-            raise ValueError("p >= 1.0")
+        if (
+            np.any(np.isnan(op))
+            or np.any(np.less_equal(op, 0.0))
+            or np.any(np.greater_equal(op, 1.0))
+        ):
+            raise ValueError("p <= 0, p >= 1 or p is NaN")
         return vec_discd_array(
             self.internal_state, irk_logseries_vec, size, op, self.lock
         )
diff --git a/mkl_random/tests/test_random.py b/mkl_random/tests/test_random.py
index 7f3b2b0..fcece59 100644
--- a/mkl_random/tests/test_random.py
+++ b/mkl_random/tests/test_random.py
@@ -926,6 +926,12 @@ def test_randomdsit_geometric(randomdist):
     np.testing.assert_array_equal(actual, desired)
 
 
+def test_geometric_exceptions():
+    for p in (0.0, -0.5, 1.5, np.nan):
+        assert_raises(ValueError, rnd.geometric, p)
+        assert_raises(ValueError, rnd.geometric, [p] * 5)
+
+
 def test_randomdist_gumbel(randomdist):
     rnd.seed(randomdist.seed, randomdist.brng)
     actual = rnd.gumbel(loc=0.123456789, scale=2.0, size=(3, 2))
@@ -1025,6 +1031,12 @@ def test_randomdist_logseries(randomdist):
     np.testing.assert_array_equal(actual, desired)
 
 
+def test_logseries_exceptions():
+    for p in (0.0, -0.5, 1.0, np.nan):
+        assert_raises(ValueError, rnd.logseries, p)
+        assert_raises(ValueError, rnd.logseries, [p] * 5)
+
+
 def test_randomdist_multinomial(randomdist):
     rs = rnd.MKLRandomState(randomdist.seed, brng=randomdist.brng)
     actual = rs.multinomial(20, [1 / 6.0] * 6, size=(3, 2))

From ed246adc41e3c92a387f24e3cd750e7bbc63319e Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Fri, 16 Oct 2026 05:25:05 -0700
Subject: [PATCH 04/10] Add gh-195 to changelog

---
 CHANGELOG.md | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index b05a04e..615ce0e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -29,8 +29,13 @@ 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)
+* `negative_binomial` now rejects `p == 0`, which MKL does not support, with a `ValueError`; the documented range of `p` is `(0, 1]` [gh-193](https://github.com/IntelPython/mkl_random/pull/193)
+* `poisson`, `negative_binomial`, `geometric` and `logseries` now reject `NaN` parameters with a `ValueError`, as NumPy does, instead of passing them to MKL [gh-193](https://github.com/IntelPython/mkl_random/pull/193)
 
 ### Fixed
+* Fixed `poisson` with a scalar `lam == 0` returning uninitialized memory as random values, since MKL requires `lam > 0`. The same applied to array-valued `lam` with `method="PTPE"`. Both now return zeros, as the default array path and NumPy already did [gh-193](https://github.com/IntelPython/mkl_random/pull/193)
+* Fixed `negative_binomial` with `p == 1` returning uninitialized memory as random values, since MKL requires `p < 1`; it now returns zeros, as NumPy does [gh-193](https://github.com/IntelPython/mkl_random/pull/193)
+* Fixed `geometric` with `p` being `NaN` returning uninitialized memory as random values [gh-193](https://github.com/IntelPython/mkl_random/pull/193)
 * 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 45e366e8009f5b36d1bbd87a09080b835eed0fc5 Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Fri, 16 Oct 2026 06:19:55 -0700
Subject: [PATCH 05/10] Make the mixed-lam check fail if MKL rejects the call

Per review: with lam = 5 for the middle element, the assertion on the zero
elements would also hold if viRngPoissonV had rejected the whole call and
left the zero-filled output buffer untouched, which is exactly the failure
mode this test is about. A lam of 1000 cannot produce a zero, so the
element now proves the call went through.

That viRngPoissonV accepts a zero lambda was checked directly, not assumed:
poisson([0, 1000, 0]) returns [0, 1005, 0] with nothing on stderr, while
the scalar and PTPE paths printed "Parameter 5 was incorrect on entry to
viRngPoisson" before this branch.
---
 mkl_random/tests/test_random.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/mkl_random/tests/test_random.py b/mkl_random/tests/test_random.py
index fcece59..30b274e 100644
--- a/mkl_random/tests/test_random.py
+++ b/mkl_random/tests/test_random.py
@@ -1326,8 +1326,11 @@ def test_poisson_zero_lam(method):
     np.testing.assert_array_equal(
         rnd.poisson(np.zeros(5), method=method), zeros
     )
-    mixed = rnd.poisson(np.array([0.0, 5.0, 0.0]), method=method)
+    # A lam too large to draw a zero, so that a call rejected as a whole,
+    # which would leave the zero-filled buffer alone, fails here
+    mixed = rnd.poisson(np.array([0.0, 1000.0, 0.0]), method=method)
     np.testing.assert_array_equal(mixed[::2], zeros[:2])
+    assert mixed[1] > 0
 
 
 def test_randomdist_power(randomdist):

From 52b8dcc39a8276a9a226bc91f365b3fc75d84ba1 Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Fri, 16 Oct 2026 06:20:00 -0700
Subject: [PATCH 06/10] Link the changelog entries to this PR

Per review: gh-193 is the unrelated set_state buffer-length issue. The
entries were written before the PR existed and kept the guessed number.
---
 CHANGELOG.md | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 615ce0e..c7654be 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -29,13 +29,13 @@ 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)
-* `negative_binomial` now rejects `p == 0`, which MKL does not support, with a `ValueError`; the documented range of `p` is `(0, 1]` [gh-193](https://github.com/IntelPython/mkl_random/pull/193)
-* `poisson`, `negative_binomial`, `geometric` and `logseries` now reject `NaN` parameters with a `ValueError`, as NumPy does, instead of passing them to MKL [gh-193](https://github.com/IntelPython/mkl_random/pull/193)
+* `negative_binomial` now rejects `p == 0`, which MKL does not support, with a `ValueError`; the documented range of `p` is `(0, 1]` [gh-195](https://github.com/IntelPython/mkl_random/pull/195)
+* `poisson`, `negative_binomial`, `geometric` and `logseries` now reject `NaN` parameters with a `ValueError`, as NumPy does, instead of passing them to MKL [gh-195](https://github.com/IntelPython/mkl_random/pull/195)
 
 ### Fixed
-* Fixed `poisson` with a scalar `lam == 0` returning uninitialized memory as random values, since MKL requires `lam > 0`. The same applied to array-valued `lam` with `method="PTPE"`. Both now return zeros, as the default array path and NumPy already did [gh-193](https://github.com/IntelPython/mkl_random/pull/193)
-* Fixed `negative_binomial` with `p == 1` returning uninitialized memory as random values, since MKL requires `p < 1`; it now returns zeros, as NumPy does [gh-193](https://github.com/IntelPython/mkl_random/pull/193)
-* Fixed `geometric` with `p` being `NaN` returning uninitialized memory as random values [gh-193](https://github.com/IntelPython/mkl_random/pull/193)
+* Fixed `poisson` with a scalar `lam == 0` returning uninitialized memory as random values, since MKL requires `lam > 0`. The same applied to array-valued `lam` with `method="PTPE"`. Both now return zeros, as the default array path and NumPy already did [gh-195](https://github.com/IntelPython/mkl_random/pull/195)
+* Fixed `negative_binomial` with `p == 1` returning uninitialized memory as random values, since MKL requires `p < 1`; it now returns zeros, as NumPy does [gh-195](https://github.com/IntelPython/mkl_random/pull/195)
+* Fixed `geometric` with `p` being `NaN` returning uninitialized memory as random values [gh-195](https://github.com/IntelPython/mkl_random/pull/195)
 * 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 02ac2f640350199854bd767072eb3d71efd4c203 Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Fri, 16 Oct 2026 06:35:34 -0700
Subject: [PATCH 07/10] Accept p == 0 in logseries, as NumPy documents

Per review: both RandomState and Generator document p in [0, 1) and return
1 there, so rejecting it was a visible difference when mkl_random is
patched into numpy.random. The sampler needs no special case, as noted in
the review: with p == 0, r = log(1 - p) is zero and the acceptance test
V >= p always holds, so irk_logseries_vec already fills ones. Checked:
logseries(0.0, 5), logseries(np.zeros(5)) and a mixed array all return ones
where p is zero, matching NumPy.

The docstring prose said 0 < p < 1 and now matches the code.
negative_binomial keeps rejecting p == 0: MKL has no such distribution
there and the number of failures has no finite answer.
---
 CHANGELOG.md                    |  1 +
 mkl_random/mklrand.pyx          | 10 +++++-----
 mkl_random/tests/test_random.py | 11 ++++++++++-
 3 files changed, 16 insertions(+), 6 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index c7654be..155597b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -36,6 +36,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 * Fixed `poisson` with a scalar `lam == 0` returning uninitialized memory as random values, since MKL requires `lam > 0`. The same applied to array-valued `lam` with `method="PTPE"`. Both now return zeros, as the default array path and NumPy already did [gh-195](https://github.com/IntelPython/mkl_random/pull/195)
 * Fixed `negative_binomial` with `p == 1` returning uninitialized memory as random values, since MKL requires `p < 1`; it now returns zeros, as NumPy does [gh-195](https://github.com/IntelPython/mkl_random/pull/195)
 * Fixed `geometric` with `p` being `NaN` returning uninitialized memory as random values [gh-195](https://github.com/IntelPython/mkl_random/pull/195)
+* Fixed `logseries` rejecting `p == 0`, which NumPy documents as valid; it now returns ones, which its sampler already produced without a special case [gh-195](https://github.com/IntelPython/mkl_random/pull/195)
 * 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)
diff --git a/mkl_random/mklrand.pyx b/mkl_random/mklrand.pyx
index 49e0106..d852d6d 100644
--- a/mkl_random/mklrand.pyx
+++ b/mkl_random/mklrand.pyx
@@ -6436,7 +6436,7 @@ cdef class _MKLRandomState:
         Draw samples from a logarithmic series distribution.
 
         Samples are drawn from a log series distribution with specified
-        shape parameter, 0 < ``p`` < 1.
+        shape parameter, 0 <= ``p`` < 1.
 
         Parameters
         ----------
@@ -6514,18 +6514,18 @@ cdef class _MKLRandomState:
         if cnp.PyArray_NDIM(op) == 0:
             fp = PyFloat_AsDouble(p)
 
-            if isnan(fp) or fp <= 0.0 or fp >= 1.0:
-                raise ValueError("p <= 0, p >= 1 or p is NaN")
+            if isnan(fp) or fp < 0.0 or fp >= 1.0:
+                raise ValueError("p < 0, p >= 1 or p is NaN")
             return vec_discd_array_sc(
                 self.internal_state, irk_logseries_vec, size, fp, self.lock
             )
 
         if (
             np.any(np.isnan(op))
-            or np.any(np.less_equal(op, 0.0))
+            or np.any(np.less(op, 0.0))
             or np.any(np.greater_equal(op, 1.0))
         ):
-            raise ValueError("p <= 0, p >= 1 or p is NaN")
+            raise ValueError("p < 0, p >= 1 or p is NaN")
         return vec_discd_array(
             self.internal_state, irk_logseries_vec, size, op, self.lock
         )
diff --git a/mkl_random/tests/test_random.py b/mkl_random/tests/test_random.py
index 30b274e..91e7279 100644
--- a/mkl_random/tests/test_random.py
+++ b/mkl_random/tests/test_random.py
@@ -1032,11 +1032,20 @@ def test_randomdist_logseries(randomdist):
 
 
 def test_logseries_exceptions():
-    for p in (0.0, -0.5, 1.0, np.nan):
+    for p in (-0.5, 1.0, np.nan):
         assert_raises(ValueError, rnd.logseries, p)
         assert_raises(ValueError, rnd.logseries, [p] * 5)
 
 
+def test_logseries_zero_p():
+    # NumPy documents p in [0, 1) and the sampler needs no special case:
+    # with p == 0 the acceptance test V >= p always holds
+    ones = np.ones(5, dtype="int32")
+    assert rnd.logseries(0.0) == 1
+    np.testing.assert_array_equal(rnd.logseries(0.0, 5), ones)
+    np.testing.assert_array_equal(rnd.logseries(np.zeros(5)), ones)
+
+
 def test_randomdist_multinomial(randomdist):
     rs = rnd.MKLRandomState(randomdist.seed, brng=randomdist.brng)
     actual = rs.multinomial(20, [1 / 6.0] * 6, size=(3, 2))

From 473d1db5f5d5eae3940ec9d024927710ade57db6 Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Fri, 16 Oct 2026 06:43:42 -0700
Subject: [PATCH 08/10] Compare dtypes in the zeros and ones checks

Same point as the review on gh-196: assert_array_equal ignores dtypes by
default, so the int32 expectations were 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 | 18 +++++++++++-------
 1 file changed, 11 insertions(+), 7 deletions(-)

diff --git a/mkl_random/tests/test_random.py b/mkl_random/tests/test_random.py
index 91e7279..cba5d45 100644
--- a/mkl_random/tests/test_random.py
+++ b/mkl_random/tests/test_random.py
@@ -1042,8 +1042,8 @@ def test_logseries_zero_p():
     # with p == 0 the acceptance test V >= p always holds
     ones = np.ones(5, dtype="int32")
     assert rnd.logseries(0.0) == 1
-    np.testing.assert_array_equal(rnd.logseries(0.0, 5), ones)
-    np.testing.assert_array_equal(rnd.logseries(np.zeros(5)), ones)
+    np.testing.assert_array_equal(rnd.logseries(0.0, 5), ones, strict=True)
+    np.testing.assert_array_equal(rnd.logseries(np.zeros(5)), ones, strict=True)
 
 
 def test_randomdist_multinomial(randomdist):
@@ -1194,9 +1194,11 @@ def test_negative_binomial_p_is_one():
     # MKL requires p < 1, while p == 1 means no trial ever fails
     zeros = np.zeros(5, dtype="int32")
     assert rnd.negative_binomial(5, 1.0) == 0
-    np.testing.assert_array_equal(rnd.negative_binomial(5, 1.0, 5), zeros)
     np.testing.assert_array_equal(
-        rnd.negative_binomial(np.full(5, 5.0), 1.0), zeros
+        rnd.negative_binomial(5, 1.0, 5), zeros, strict=True
+    )
+    np.testing.assert_array_equal(
+        rnd.negative_binomial(np.full(5, 5.0), 1.0), zeros, strict=True
     )
 
 
@@ -1331,14 +1333,16 @@ def test_poisson_zero_lam(method):
     # MKL requires lam > 0, while lam == 0 always gives zero events
     zeros = np.zeros(5, dtype="int32")
     assert rnd.poisson(0.0, method=method) == 0
-    np.testing.assert_array_equal(rnd.poisson(0.0, 5, method=method), zeros)
     np.testing.assert_array_equal(
-        rnd.poisson(np.zeros(5), method=method), zeros
+        rnd.poisson(0.0, 5, method=method), zeros, strict=True
+    )
+    np.testing.assert_array_equal(
+        rnd.poisson(np.zeros(5), method=method), zeros, strict=True
     )
     # A lam too large to draw a zero, so that a call rejected as a whole,
     # which would leave the zero-filled buffer alone, fails here
     mixed = rnd.poisson(np.array([0.0, 1000.0, 0.0]), method=method)
-    np.testing.assert_array_equal(mixed[::2], zeros[:2])
+    np.testing.assert_array_equal(mixed[::2], zeros[:2], strict=True)
     assert mixed[1] > 0
 
 

From f1692a2a6320d4cd52d31978165941033480e2a6 Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Fri, 16 Oct 2026 09:28:51 -0700
Subject: [PATCH 09/10] Say uninitialized, not zero-filled, in the poisson test

Per review: vec_Poisson_array allocates the output with PyArray_SimpleNew,
and the scalar path with np.empty, so a rejected call leaves whatever was
in memory rather than zeros. The check itself is unaffected.
---
 mkl_random/tests/test_random.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/mkl_random/tests/test_random.py b/mkl_random/tests/test_random.py
index cba5d45..be50a25 100644
--- a/mkl_random/tests/test_random.py
+++ b/mkl_random/tests/test_random.py
@@ -1340,7 +1340,7 @@ def test_poisson_zero_lam(method):
         rnd.poisson(np.zeros(5), method=method), zeros, strict=True
     )
     # A lam too large to draw a zero, so that a call rejected as a whole,
-    # which would leave the zero-filled buffer alone, fails here
+    # which would leave the uninitialized buffer alone, fails here
     mixed = rnd.poisson(np.array([0.0, 1000.0, 0.0]), method=method)
     np.testing.assert_array_equal(mixed[::2], zeros[:2], strict=True)
     assert mixed[1] > 0

From 0279e2c1526f582a32a37edee692044b9e4182fc Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Fri, 16 Oct 2026 09:29:12 -0700
Subject: [PATCH 10/10] Document p in the logseries docstring

Per review, older than this PR but in the lines being touched: the
Parameters block described loc and scale, which the function does not take,
and never mentioned p, while Returns claimed the values lie in [0, n]. The
support is the positive integers.

gh-165 rewrites this block as well, so that PR keeps the final wording; the
range there needs to read [0, 1) after this one.
---
 mkl_random/mklrand.pyx | 8 +++-----
 1 file changed, 3 insertions(+), 5 deletions(-)

diff --git a/mkl_random/mklrand.pyx b/mkl_random/mklrand.pyx
index d852d6d..f8d09c1 100644
--- a/mkl_random/mklrand.pyx
+++ b/mkl_random/mklrand.pyx
@@ -6440,10 +6440,8 @@ cdef class _MKLRandomState:
 
         Parameters
         ----------
-        loc : float
-
-        scale : float > 0.
-
+        p : float or array_like of floats
+            Shape parameter, 0 <= ``p`` < 1.
         size : int or tuple of ints, optional
             Output shape.  If the given shape is, e.g., ``(m, n, k)``, then
             ``m * n * k`` samples are drawn.  Default is None, in which case a
@@ -6452,7 +6450,7 @@ cdef class _MKLRandomState:
         Returns
         -------
         samples : ndarray or scalar
-                  where the values are all integers in  [0, n].
+                  where the values are all positive integers.
 
         See Also
         --------