Add discrete distributions to `dpnp.random.Generator` by ndgrigorian · Pull Request #3159 · IntelPython/dpnp

ndgrigorian · 2026-11-04T19:46:03Z

This PR adds integers and the discrete distributions to dpnp.random.Generator, as part of the series splitting #3147. It is stacked on #3158

Adds

integers, using NumPy's bounded integer algorithm (Lemire's method)
poisson, binomial (inversion and BTPE), negative_binomial, geometric, logseries, zipf and hypergeometric (HRUA)
noncentral_chisquare and noncentral_f, which need poisson

The discrete samplers return int64 and compute in double precision, or in single precision on devices without fp64. In single precision, host parameters are clipped so that they stay inside the domain once rounded, and the large log-factorial terms of BTPE and HRUA are gathered into one log1p so that they don't cancel. zipf uses expm1/log1p, deliberately departing from NumPy, whose tail for a close to 1 is too heavy above 2**53 (documented in random.rst)

As in #3158, every rejection sampler returns the smallest int64 for device parameters that would make it loop forever

Also adds umul_hi, which replaces sycl::mul_hi; DPC++ 2026.2 miscompiles sycl::mul_hi on Intel GPUs when one operand is a compile-time constant

Have you provided a meaningful PR description?
Have you added a test, reproducer or referred to an issue with a reproducer?
Have you tested your changes locally for CPU and GPU devices?
Have you made sure that new changes do not introduce compiler warnings?
Have you checked performance impact of proposed changes?
Have you added documentation for your changes, if necessary?
Have you added your changes to the changelog?

github-actions · 2026-11-04T20:30:35Z

View rendered docs @ https://intelpython.github.io/dpnp/pull/3159/index.html

github-actions · 2026-11-04T20:39:49Z

Array API standard conformance tests for dpnp=0.22.0dev1=np2py314h8d9cdd5_20 ran successfully.
Passed: 1491
Failed: 0
Skipped: 6

coveralls · 2026-11-04T21:36:58Z

Coverage is 84.705% — feature/device-rng-discrete into feature/device-rng-continuous. No base build found for feature/device-rng-continuous.

Add integers and the poisson, binomial, negative_binomial, geometric, logseries, zipf and hypergeometric samplers, along with noncentral_chisquare and noncentral_f, which need poisson. The discrete samplers return int64 and compute in double precision, or in single precision on devices without fp64, where host parameters are clipped to stay inside their domain once rounded.

vlad-perevezentsev · 2026-11-12T18:42:31Z

+std::int64_t binomial_inversion(BitGenT &bg, const std::int64_t n, const T p)
+{
+    const T q = T(1) - p;
+    const T qn = sycl::exp(static_cast<T>(n) * sycl::log(q));


Blocking: On the float32 path (devices without fp64), q = 1 - p is rounded before the log, so qn is wrong for small p. All draws are 0 for p below about 3e-8, and the upper tail is cut when 1 - p rounds up:

import numpy, dpnp from dpnp.random import Generator, default_rng # the single-precision path of devices without fp64 Generator._discrete_float_type = lambda self: numpy.dtype(numpy.float32) rng = default_rng(0) print(rng.binomial(3 * 10**9, 1e-8, size=10**6).mean()) # 0.0, expected 30 print(rng.binomial(21428571, 1.4e-6, size=10**6).mean()) # 25.94, expected 30

sycl::exp(n * sycl::log1p(-p)) fixes both in a float32 emulation of this loop. Could you also add such a case to test_pmf?

vlad-perevezentsev · 2026-11-12T18:42:32Z

+    return request.param
+
+    @pytest.mark.parametrize("usm_type", ["shared", "host"])
+    def test_usm_type_follows_bounds(self, usm_type):


test_usm_type_follows_bounds is indented into the discrete_precision fixture, after its return, and test_device_integer_param_dtype (line 1417) into _assert_same_error, so pytest collects neither. Moved into TestIntegers and TestDiscreteMethods, the first one passes, but the second one fails: rng.binomial(n * 1j, 0.5) returns zeros instead of raising TypeError (see the comment on _as_int64).

vlad-perevezentsev · 2026-11-12T18:42:32Z

+def _as_int64(x):
+    """Truncate an integer parameter to int64, as NumPy does."""
+
+    if dpnp.is_supported_array_type(x):
+        return dpnp.asarray(x).astype(_int64)
+    if numpy.ndim(x) > 0:
+        return numpy.asarray(x).astype(_int64)


_as_int64 casts before any check. So a complex device array is cast silently (the integer branch of _check_device_cast only ever sees int64), and NumPy float or uint64 arrays are truncated, where NumPy raises TypeError:

import numpy, dpnp rng = dpnp.random.default_rng(0) print(rng.binomial(dpnp.full(4, 5.0) * 1j, 0.5)) # [0 0 0 0], no TypeError print(rng.binomial(numpy.array([5.5]), 0.5)) # [1]; NumPy raises TypeError

Could it call _check_device_cast(x, _int64) before astype, and check numpy.can_cast(x.dtype, numpy.int64) for NumPy arrays?

vlad-perevezentsev · 2026-11-12T18:42:32Z

+    }
+
+    // inversion
+    const T z = sycl::ceil(-standard_exponential<T>(bg) / sycl::log1p(-p));


On the float32 path, geometric returns 0, which is outside its support, when the uniform draw is exactly 0 (ceil(0) == 0). That is once in 2**24 draws:

import numpy, dpnp from dpnp.random import Generator, default_rng # the single-precision path of devices without fp64 Generator._discrete_float_type = lambda self: numpy.dtype(numpy.float32) rng = default_rng(1) print(sum(int((rng.geometric(0.1, size=2**25) == 0).sum()) for _ in range(8))) # 13

Could it take sycl::fmax(z, T(1))?

vlad-perevezentsev · 2026-11-12T18:42:32Z

+
+        low = int(low)
+        high = int(high)
+        _check_bounds(low, high)


NumPy returns an empty array, without checking the bounds, when size has no elements. Here integers raises:

import numpy, dpnp print(numpy.random.default_rng(0).integers(5, 5, size=0)) # [] dpnp.random.default_rng(0).integers(5, 5, size=0) # ValueError: low >= high

Could _integers return an empty array first in that case, for the array bounds too?

vlad-perevezentsev · 2026-11-12T18:42:32Z

+    def poisson(self, lam=1.0, size=None, *, usm_type=None):
+        """
+        Draw samples from a Poisson distribution.
+
+        For full documentation refer to :obj:`numpy.random.Generator.poisson`.
+
+        Limitations
+        -----------
+        `usm_type` is not in NumPy's signature.
+
+        """
+
+        self._check_queues(lam)
+        _check_constraint(


The new methods (poisson to integers, and the two noncentral ones) have no Parameters and Returns sections, which dpnp/AGENTS.md asks for. So the rendered docs say nothing about usm_type or the int64 result. random and normal show the layout.

vlad-perevezentsev · 2026-11-12T18:42:32Z

+    if (!(a > T(1))) {
+        // invalid, and every draw would be rejected
+        return 0;


Nit: For an invalid a on the device, zipf returns 0, while the other discrete samplers return the smallest int64, as the new line in dpnp/libtensor/AGENTS.md says:

import dpnp rng = dpnp.random.default_rng(0) print(rng.zipf(dpnp.asarray([0.5, dpnp.nan]))) # [0 0] print(rng.geometric(dpnp.asarray([dpnp.nan]))) # [-9223372036854775808]

Could it return std::numeric_limits<std::int64_t>::min() too, with a line in test_device_params_terminate?

vlad-perevezentsev · 2026-11-12T18:42:32Z

+            def _to_bound(x, shift):
+                if dpnp.is_supported_array_type(x):
+                    x = dpnp.asarray(x, dtype=wide, sycl_queue=queue)
+                    return (x - shift if shift else x).astype(dtype)


Nit: astype copies by default, so an int64 bound on the device is copied here once more after x - shift, and _as_int64 (line 306) copies an int64 device array as well. astype(..., copy=False) in both places avoids these device copies.

vlad-perevezentsev · 2026-11-12T18:42:32Z

+// sycl::mul_hi is not used: optimized DPC++ 2026.2 builds get it wrong on
+// Intel GPUs when one 64-bit operand is a compile-time constant (all results
+// wrong at -O2, right at -O0). oneMKL's device RNG uses 32-bit halves too.


Nit: This works around a DPC++ bug, and AGENTS.md asks for the full issue URL in such a comment, so it can be reverted later. Is there an intel/llvm issue for the sycl::mul_hi miscompile? The new line in dpnp/libtensor/AGENTS.md could link it too.

vlad-perevezentsev · 2026-11-12T18:42:32Z

+        else -math.inf
+    ),
+    "zipf": lambda k, a: (
+        -a * math.log(k) - math.log(_zeta(a)) if k >= 1 else -math.inf


Nit: The zipf pmf calls _zeta(a), a sum of 10**5 terms, for every k, so the four zipf test_pmf cases take about 9.5 s on the CPU, most of TestDiscreteMethods. @functools.lru_cache on _zeta would remove that.

ndgrigorian requested review from antonwolfy and vlad-perevezentsev as code owners November 4, 2026 19:46

This was referenced Nov 4, 2026

Add dpnp.random.Generator on the oneMKL device RNG API #3157

Open

Add multivariate distributions and permutations to dpnp.random.Generator #3160

Open

ndgrigorian force-pushed the feature/device-rng-discrete branch from 9fa550e to bb0eae2 Compare November 4, 2026 19:47

ndgrigorian force-pushed the feature/device-rng-continuous branch 2 times, most recently from 0160558 to adf7178 Compare November 4, 2026 19:50

ndgrigorian force-pushed the feature/device-rng-discrete branch from bb0eae2 to 7918e32 Compare November 4, 2026 19:50

ndgrigorian force-pushed the feature/device-rng-continuous branch from adf7178 to cfa00bd Compare November 6, 2026 17:55

ndgrigorian force-pushed the feature/device-rng-discrete branch from 7918e32 to 66e7a74 Compare November 6, 2026 17:55

ndgrigorian added 2 commits November 8, 2026 13:03

Pass the dpnp array straight to assert_allclose in the poisson test

f48c80e

ndgrigorian force-pushed the feature/device-rng-discrete branch from 6388643 to f48c80e Compare November 8, 2026 21:54

ndgrigorian force-pushed the feature/device-rng-continuous branch from df5cb09 to 7c8127b Compare November 8, 2026 21:54

vlad-perevezentsev reviewed Nov 12, 2026

View reviewed changes

Conversation

ndgrigorian commented Nov 4, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

github-actions Bot commented Nov 4, 2026

Uh oh!

github-actions Bot commented Nov 4, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

coveralls commented Nov 4, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

This branch has not been deployed

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

ndgrigorian commented Nov 4, 2026 •

edited

Loading

github-actions Bot commented Nov 4, 2026 •

edited

Loading

coveralls commented Nov 4, 2026 •

edited

Loading