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

CHANGELOG.md

@@ -12,6 +12,7 @@ This release drops support for Python 3.10, making Python 3.11 the minimum requi
 
 * Added support for boolean input in `dpnp.positive` and the `+` operator of `dpnp.ndarray` [#3128](https://github.com/IntelPython/dpnp/pull/3128)
 * Added `dpnp.top_k` to return the `k` largest or smallest elements and their indices along an axis, with a `sorted` keyword to return them in sorted order [#3133](https://github.com/IntelPython/dpnp/pull/3133)
+* Added the `dpnp.random.Generator` API with `default_rng`, `BitGenerator`, `SeedSequence` and the `Philox4x3210` bit generator, sampling on the device with the device RNG API of oneMKL or oneMath [#3157](https://github.com/IntelPython/dpnp/pull/3157)
 
 ### Changed
 

doc/known_words.txt

@@ -3,6 +3,7 @@ ary
 backend
 bandlimited
 bincount
+BitGenerator
 bitwise
 Blackman
 boolean
@@ -79,19 +80,24 @@ ord
 orthonormal
 radix
 Penrose
+Philox
 Polyutils
 pre
 prepend
 prepended
 prepending
+reentrant
 representable
 resampling
+rng
 runtimes
 scikit
 se
+SeedSequence
 signbit
 signum
 sinc
+strided
 subarray
 subarrays
 subclasses
@@ -118,4 +124,5 @@ Vetterline
 von
 Weibull
 whitespace
+ziggurat
 Zipf

doc/reference/random.rst

@@ -3,10 +3,121 @@
 Random sampling (:mod:`dpnp.random`)
 ====================================
 
+.. hint:: `NumPy API Reference: Random sampling <https://numpy.org/doc/stable/reference/random/index.html>`_
+
+DPNP draws random numbers on the SYCL device, with the counter-based
+Philox4x32x10 generator of the oneMKL device RNG API (oneMath's in a
+oneMath build). It provides the :class:`Generator` API of NumPy,
+recommended for new code. The legacy :class:`RandomState` API and the
+module-level functions still use the host-side oneMKL RNG.
+
+.. code-block:: python
+
+    import dpnp
+
+    rng = dpnp.random.default_rng(seed=42)
+    x = rng.standard_normal((3, 4))                # on the default device
+    y = rng.random(5, usm_type="shared")
+
+    rng = dpnp.random.default_rng(42, device="cpu")
+
+A generator draws all its samples on one SYCL queue, chosen with the
+``device`` or ``sycl_queue`` keyword when it is created. All samples are
+:class:`dpnp.ndarray`, zero-dimensional when ``size`` is ``None``.
+
+Differences from NumPy
+----------------------
+
+* The streams differ from NumPy's. For a given seed, output size and data
+  type, the raw bits and the uniform floats are the same on every device,
+  while the normal samples agree to a few ULP.
+* The values depend on the output size: a draw is not a prefix of a
+  larger one, and two draws of ``n`` values differ from one draw of
+  ``2 * n``. Code splitting its draws into chunks gets other values.
+* The floating-point distributions take a ``dtype`` keyword, as in CuPy,
+  which defaults to the default floating-point type of the device, so on a
+  device without double precision the same seed gives other values by
+  default.
+* All methods take a ``usm_type`` keyword. With ``usm_type=None``, the
+  default, the samples take the USM type of the array parameters, or
+  ``"device"`` when there are none.
+* Parameters on the host (Python scalars, lists and NumPy arrays) are
+  checked as in NumPy. Parameters in :class:`dpnp.ndarray` are
+  not, since that would wait for the device.
+  Invalid ones give undefined values, such as NaN.
+* Parameters given as lists or NumPy arrays are copied to the device first,
+  which waits for the work already submitted to the queue. Scalars and
+  :class:`dpnp.ndarray` parameters don't wait.
+* Generators and bit generators cannot be pickled or deep-copied, and the
+  bit generators have no ``advance`` or ``jumped`` method.
+
+Random Generator
+----------------
+
+.. autosummary::
+   :toctree: generated/
+   :nosignatures:
+
+   default_rng
+   Generator
+
+Accessing the bit generator and spawning
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. autosummary::
+   :toctree: generated/
+   :nosignatures:
+
+   Generator.bit_generator
+   Generator.spawn
+   Generator.sycl_device
+   Generator.sycl_queue
+
+Simple random data
+^^^^^^^^^^^^^^^^^^
+
+.. autosummary::
+   :toctree: generated/
+   :nosignatures:
+
+   Generator.random
+
+Distributions
+^^^^^^^^^^^^^
+
+.. autosummary::
+   :toctree: generated/
+   :nosignatures:
+
+   Generator.normal
+   Generator.standard_normal
+   Generator.uniform
+
+Bit generators
+--------------
+
+.. autosummary::
+   :toctree: generated/
+   :nosignatures:
+
+   BitGenerator
+   BitGenerator.lock
+   BitGenerator.random_raw
+   BitGenerator.seed_seq
+   BitGenerator.spawn
+   BitGenerator.state
+   BitGenerator.sycl_device
+   BitGenerator.sycl_queue
+   Philox4x3210
+   SeedSequence
+
+Legacy random generation
+------------------------
+
 .. hint:: `NumPy API Reference: Legacy random generation <https://numpy.org/doc/stable/reference/random/legacy.html>`_
 
 Simple random data
-------------------
+^^^^^^^^^^^^^^^^^^
 
 .. autosummary::
    :toctree: generated/
@@ -23,9 +134,8 @@ Simple random data
    choice
    bytes
 
-
 Permutations
-------------
+^^^^^^^^^^^^
 
 .. autosummary::
    :toctree: generated/
@@ -34,9 +144,8 @@ Permutations
    shuffle
    permutation
 
-
 Distributions
--------------
+^^^^^^^^^^^^^
 
 .. autosummary::
    :toctree: generated/
@@ -78,9 +187,8 @@ Distributions
    weibull
    zipf
 
-
-Random generator
-----------------
+RandomState
+^^^^^^^^^^^
 
 .. autosummary::
    :toctree: generated/

dpnp/AGENTS.md

@@ -6,6 +6,7 @@
 - `dpnp_algo/dpnp_elementwise_common.py`: ufunc wrapper classes (`DPNPUnaryFunc`, `DPNPBinaryFunc`, two-output and specialised variants)
 - `dpnp_utils/`: shared helpers (`get_usm_allocations`, `map_dtype_to_device`, reduction wrappers)
 - `linalg/`, `fft/`, `random/`, `scipy/`, `lib/`: subpackages, each exporting through its own `__init__.py`
+- `random/`: `_bit_generator.py` (`BitGenerator`, `Philox4x3210`; `SeedSequence` is NumPy's) and `_generator.py` (`Generator`, which implements every distribution). The kernels are in `dpnp.libtensor._rng_impl`.
 - Legacy, don't extend: the Cython `dpnp_algo` path and `dpnp_backend_c` (still used by `dpnp.partition` and `dpnp.random`), and the NumPy fallback `call_origin`. New kernels go into `dpnp/libtensor` (see `dpnp/libtensor/AGENTS.md`).
 
 ## Adding a Public Function
@@ -20,6 +21,7 @@
 - Extension calls are ordered with `_manager = dpctl.utils.SequentialOrderManager[exec_q]`: pass `depends=_manager.submitted_events`, then `_manager.add_event_pair(ht_ev, ev)`. Don't add blocking waits.
 - Use `dpnp.default_float_type(...)` or `map_dtype_to_device(...)` for the default float type.
 - Deprecated ufuncs use `DPNPDeprecatedUnaryFunc(..., deprecated_msg=...)`.
+- A `dpnp.random` distribution is a `Generator` method: it calls `self._check_queues(...)` on its array parameters (raises `ExecutionPlacementError`), checks host parameters as NumPy does, and draws through `_cont` or `_float_sample` into `BitGenerator._fill`. Parameters on the device are not checked, as the check waits for the device. With `usm_type=None`, the output and the temporaries take the coerced USM type of the array parameters (`Generator._sample`).
 - Unsupported NumPy keywords raise `NotImplementedError` (`dpnp.check_limitations`) and are listed under `Limitations` in the docstring.
 - The Array API spec version is set in three places that change together: `__array_api_version__` in `dpnp/tensor/_array_api.py`, `ARRAY_API_TESTS_VERSION` in `.github/workflows/conda-package.yml` and `doc/reference/array_api.rst` (plus the version test in `dpnp/tests/test_ndarray.py`). Inspection utilities come from `__array_namespace_info__` in `dpnp_array_api_info.py`.
 - `dpnp.tensor` replaces `dpctl.tensor`, which no longer exists here; `SyclQueue`, `SyclDevice`, `dpctl.memory` and `SequentialOrderManager` still come from `dpctl`.

dpnp/CMakeLists.txt

@@ -229,6 +229,7 @@ add_subdirectory(libtensor/source)
 add_subdirectory(libtensor/extensions/blas)
 add_subdirectory(libtensor/extensions/fft)
 add_subdirectory(libtensor/extensions/lapack)
+add_subdirectory(libtensor/extensions/rng)
 add_subdirectory(libtensor/extensions/vm)
 add_subdirectory(libtensor/extensions/sparse)
 

dpnp/libtensor/AGENTS.md

@@ -3,7 +3,7 @@
 ## Layout
 - `include/kernels/`: header-only SYCL kernels in namespaces under `dpnp::tensor::kernels`; `include/utils/`: type dispatch, type, offset and memory-overlap helpers (`type_dispatch.hpp`, `type_dispatch_building.hpp`, `type_utils.hpp`, `memory_overlap.hpp`, ...); `include/ext/`: dispatch helpers (`init_dispatch_vector`, `init_dispatch_table`, `DispatchTable`) and validation helpers (`validation_utils.hpp`). The headers ship with dpnp; `dpnp.get_include()` and `python -m dpnp --tensor-include-dir` return `include/`.
 - `source/`: bindings built by `source/CMakeLists.txt` into `_tensor_impl`, `_tensor_elementwise_impl`, `_tensor_reductions_impl`, `_tensor_sorting_impl`, `_tensor_accumulation_impl`, `_tensor_linalg_impl`, `_statistics_impl` and `_window_impl`. The functions dpnp adds beyond the Array API are listed in the `_dpnp_*_sources` variables there and compiled into the same modules.
-- `extensions/{blas,fft,lapack,sparse,vm}/`: oneMKL/oneMath extensions, one module `_<name>_impl` per directory.
+- `extensions/{blas,fft,lapack,rng,sparse,vm}/`: oneMKL/oneMath extensions, one module `_<name>_impl` per directory.
 - Every module installs into `dpnp/libtensor` and is imported as `dpnp.libtensor._<name>_impl`. A new kernel goes into `source/` and an existing module; only code that needs oneMKL gets a new extension.
 
 ## Adding a ufunc Kernel (model: `fabs`)
@@ -21,8 +21,9 @@ A new extension's `CMakeLists.txt` lists its sources and calls `build_dpnp_pybin
 - Validate inputs before launching: queue compatibility, contiguity, writability, memory overlap (including the same logical tensor), dtype and ndim. Reuse the helpers neighbouring files use: `check_queue`, `check_c_contig`, `check_no_overlap`, `check_writable` from `include/ext/validation_utils.hpp` (`putmask`, `statistics`, `isclose`, `interpolate`), or `dpnp::utils::queues_are_compatible`, `dpnp::tensor::overlap::MemoryOverlap` and `dpnp::tensor::validation::CheckWritable` (extensions and the Array API functions). Don't silently accept strided input.
 - Throw `py::value_error` for bad arguments and `std::runtime_error` for internal failures; catch oneMKL and SYCL exceptions and rethrow them.
 - Guard double, complex<double> and half with `dpnp::tensor::type_utils::validate_type_for_device<T>(q)`; don't use double literals or casts in kernels that must run without fp64.
-- `blas`, `fft`, `lapack` and `sparse` build against either oneMKL or oneMath (`USE_ONEMATH`); `vm` has kernels only with oneMKL: with oneMath it builds empty and `_is_available()` returns `False`.
+- `blas`, `fft`, `lapack` and `sparse` build against either oneMKL or oneMath (`USE_ONEMATH`); `vm` has kernels only with oneMKL: with oneMath it builds empty and `_is_available()` returns `False`; `rng` uses the header-only device RNG API of oneMKL, or of oneMath with `USE_ONEMATH`, and links neither library.
 - Release the GIL around blocking oneMKL calls (`py::gil_scoped_release`).
+- `rng` draws depend only on the engine state and the number of elements, never on the device or the layout chosen for it (`extensions/rng/kernels/sample.hpp`): raw bits and uniform floats are bit-identical across devices, samplers using transcendental functions agree to a few ULP. Each engine's `items_per_work_item` (`kernels/engines.hpp`) is part of its stream definition.
 - An extension whose kernels may be launched over more than `INT_MAX` work-items passes `LARGE_RANGES` to `build_dpnp_pybind11_ext`, which adds `-fno-sycl-id-queries-fit-in-int` at compile and link time (DPC++ 2026.2 and later need it at link time). `source/CMakeLists.txt` and `dpnp_backend_c` pass it to both `target_compile_options` and `target_link_options`.
 - clang-format 23 formats C++ (`.clang-format`: LLVM base, 4-space indent, 80 columns); gersemi formats CMake at 88 columns, except `dpnp/backend/cmake/Modules/`.
 

dpnp/libtensor/extensions/rng/CMakeLists.txt

@@ -0,0 +1,66 @@
+# -*- coding: utf-8 -*-
+# *****************************************************************************
+# Copyright (c) 2026, Intel Corporation
+# All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+# - Redistributions of source code must retain the above copyright notice,
+#   this list of conditions and the following disclaimer.
+# - Redistributions in binary form must reproduce the above copyright notice,
+#   this list of conditions and the following disclaimer in the documentation
+#   and/or other materials provided with the distribution.
+# - Neither the name of the copyright holder nor the names of its contributors
+#   may be used to endorse or promote products derived from this software
+#   without specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
+# THE POSSIBILITY OF SUCH DAMAGE.
+# *****************************************************************************
+
+set(python_module_name _rng_impl)
+set(_module_src
+    ${CMAKE_CURRENT_SOURCE_DIR}/normal.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/rng_py.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/sample.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/uniform.cpp
+)
+
+# a launch takes one work-item per 8 elements, which can exceed INT_MAX
+build_dpnp_pybind11_ext(
+    ${python_module_name}
+    SOURCES ${_module_src}
+    LARGE_RANGES
+)
+
+# the samplers rely on IEEE rounding, which reassociation would change
+if(WIN32)
+    target_compile_options(${python_module_name} PRIVATE /clang:-fno-fast-math)
+else()
+    target_compile_options(${python_module_name} PRIVATE -fno-fast-math)
+endif()
+
+# the device RNG API is header-only, so no oneMKL or oneMath library is linked
+if(_use_onemath)
+    target_compile_options(${python_module_name} PRIVATE -DUSE_ONEMATH)
+    target_include_directories(
+        ${python_module_name}
+        SYSTEM
+        PRIVATE $<TARGET_PROPERTY:${ONEMATH_LIB},INTERFACE_INCLUDE_DIRECTORIES>
+    )
+else()
+    target_include_directories(
+        ${python_module_name}
+        SYSTEM
+        PRIVATE $<TARGET_PROPERTY:MKL::MKL_SYCL,INTERFACE_INCLUDE_DIRECTORIES>
+    )
+endif()