Reimplement dpnp.random using oneMKL device RNG API

CHANGELOG.md

@@ -10,6 +10,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 * 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, on the oneMKL device RNG API [#3147](https://github.com/IntelPython/dpnp/pull/3147)
+* Added the `dpnp.random.PCG64DXSM` bit generator, which matches `numpy.random.PCG64DXSM` for raw values, `random`, `uniform` and `bytes` [#3147](https://github.com/IntelPython/dpnp/pull/3147)
+* Added `dpnp.random.get_bit_generator` and `dpnp.random.set_bit_generator`, and support for a bit generator as the seed of `dpnp.random.RandomState` [#3147](https://github.com/IntelPython/dpnp/pull/3147)
+* Added the `dpnp.random.MRG32k3a` bit generator [#3147](https://github.com/IntelPython/dpnp/pull/3147)
 
 ### Changed
 
@@ -18,6 +22,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 * Merged the `ufunc` and `indexing` extensions into `dpnp.libtensor._tensor_elementwise_impl` and `dpnp.libtensor._tensor_impl`, while `statistics` and `window` extensions became part of libtensor [#3131](https://github.com/IntelPython/dpnp/pull/3131)
 * Moved `blas`, `fft`, `lapack`, `sparse` and `vm` extensions into `dpnp/libtensor/extensions` [#3131](https://github.com/IntelPython/dpnp/pull/3131)
 * Moved the libtensor headers shipped with dpnp from `dpnp/tensor/libtensor/include` to `dpnp/libtensor/include`, and changed `dpnp.get_include` to return that directory [#3131](https://github.com/IntelPython/dpnp/pull/3131)
+* Reimplemented `dpnp.random.RandomState` and the `dpnp.random` functions on device kernels, giving the same values on every device. `device`, `usm_type` and `sycl_queue` are now keyword-only, `randint(dtype=int)` returns `int64` and `shuffle` requires a `dpnp.ndarray` [#3147](https://github.com/IntelPython/dpnp/pull/3147)
 
 ### Deprecated
 
@@ -29,6 +34,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Removed
 
 * Removed the unused `FourZeroOffsets_Indexer`, `FixedDimStridedIndexer` and `ThreeOffsets_FixedDimStridedIndexer` from the libtensor header `utils/offset_utils.hpp` [#3131](https://github.com/IntelPython/dpnp/pull/3131)
+* Removed the legacy backend implementation of `dpnp.random`, which no longer needs `onemkl-sycl-rng` [#3147](https://github.com/IntelPython/dpnp/pull/3147)
+* Removed the unused parts of the legacy `dpnp_backend_c` library, including its empty GTest suite, the `DPNP_BACKEND_TESTS` CMake option and the NumPy fallback helper `call_origin` [#3147](https://github.com/IntelPython/dpnp/pull/3147)
 
 ### Fixed
 

CMakeLists.txt

@@ -42,7 +42,6 @@ option(
     "Build dpnp tensor pybind11 offloading extensions with coverage instrumentation"
     OFF
 )
-option(DPNP_BACKEND_TESTS "Enable building of DPNP backend test suite" OFF)
 option(
     DPNP_WITH_REDIST
     "Build DPNP assuming DPC++ redistributable is installed into Python prefix"

conda-recipe-cf/meta.yaml

@@ -56,7 +56,6 @@ requirements:
     - {{ pin_compatible('onemkl-sycl-blas', min_pin='x.x', max_pin='x') }}
     - {{ pin_compatible('onemkl-sycl-dft', min_pin='x.x', max_pin='x') }}
     - {{ pin_compatible('onemkl-sycl-lapack', min_pin='x.x', max_pin='x') }}
-    - {{ pin_compatible('onemkl-sycl-rng', min_pin='x.x', max_pin='x') }}
     - {{ pin_compatible('onemkl-sycl-vm', min_pin='x.x', max_pin='x') }}
     - numpy
     - ocl-icd-system  # [linux]

conda-recipe/meta.yaml

@@ -69,7 +69,6 @@ requirements:
       - {{ pin_compatible('onemkl-sycl-blas', min_pin='x.x', max_pin='x') }}
       - {{ pin_compatible('onemkl-sycl-dft', min_pin='x.x', max_pin='x') }}
       - {{ pin_compatible('onemkl-sycl-lapack', min_pin='x.x', max_pin='x') }}
-      - {{ pin_compatible('onemkl-sycl-rng', min_pin='x.x', max_pin='x') }}
       - {{ pin_compatible('onemkl-sycl-vm', min_pin='x.x', max_pin='x') }}
       - {{ pin_compatible('onemkl-sycl-sparse', min_pin='x.x', max_pin='x') }}
       - numpy

conda-recipe/rattler_recipe.yaml

@@ -70,7 +70,6 @@ requirements:
     - ${{ pin_compatible('onemkl-sycl-blas', lower_bound='x.x', upper_bound='x') }}
     - ${{ pin_compatible('onemkl-sycl-dft', lower_bound='x.x', upper_bound='x') }}
     - ${{ pin_compatible('onemkl-sycl-lapack', lower_bound='x.x', upper_bound='x') }}
-    - ${{ pin_compatible('onemkl-sycl-rng', lower_bound='x.x', upper_bound='x') }}
     - ${{ pin_compatible('onemkl-sycl-vm', lower_bound='x.x', upper_bound='x') }}
     - ${{ pin_compatible('onemkl-sycl-sparse', lower_bound='x.x', upper_bound='x') }}
     - numpy

doc/known_words.txt

@@ -42,9 +42,11 @@ Hanning
 histogrammed
 Hypergeometric
 kwargs
+Kelton
 iaxis
 iinfo
 Infs
+L'Ecuyer
 intp
 ints
 iterable
@@ -86,15 +88,19 @@ prepended
 prepending
 representable
 resampling
+RngStreams
 runtimes
 scikit
 se
 signbit
 signum
+Simard
 sinc
 subarray
 subarrays
 subclasses
+substream
+substreams
 subtype
 SyclDevice
 SyclQueue

doc/reference/random.rst

@@ -3,10 +3,185 @@
 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, by default with the
+counter-based Philox4x32x10 generator of the oneMKL device RNG API, or with
+NumPy's PCG64DXSM generator or the MRG32k3a generator of the oneMKL device
+RNG API. It provides both NumPy random APIs: the :class:`Generator` API recommended for new code,
+and the legacy :class:`RandomState` API with the module-level functions
+built on it.
+
+.. code-block:: python
+
+    import dpnp
+
+    rng = dpnp.random.default_rng(seed=42)
+    x = rng.standard_normal((3, 4))                # on the default device
+    y = rng.integers(0, 10, size=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, and the samples
+depend only on the seed, not on the device, except that the discrete
+distributions also depend on its support of double precision. All samples are
+:class:`dpnp.ndarray`, zero-dimensional when ``size`` is ``None``.
+
+Differences from NumPy
+----------------------
+
+* The same seed gives the same stream on every device, but for the discrete
+  distributions noted below. The streams differ from NumPy's, except that
+  :class:`PCG64DXSM` gives NumPy's raw values, :meth:`Generator.random`,
+  :meth:`Generator.uniform` and :meth:`Generator.bytes`, and can exchange
+  states with :class:`numpy.random.PCG64DXSM`.
+* The floating-point distributions take a ``dtype`` keyword, as in CuPy,
+  which defaults to the default floating-point type of the device. All
+  methods take a ``usm_type`` keyword.
+* Parameters on the host (Python scalars, lists and NumPy arrays) are
+  checked as in NumPy. In :class:`Generator`, parameters in
+  :class:`dpnp.ndarray` are not, since that would wait for the device.
+  Invalid ones give undefined values, such as NaN or the smallest ``int64``,
+  but the kernels handle them safely: they never loop without end or access
+  memory out of bounds. Floating-point arrays given for integer parameters
+  are converted as by :meth:`dpnp.ndarray.astype`, whose result for NaN and
+  out-of-range values is undefined.
+* The discrete distributions return ``int64`` values. They compute in double
+  precision on devices supporting it, and in single precision otherwise.
+  There, the streams differ from the ones of devices with double precision,
+  and some samples are only multiples of a power of two: the ones above
+  ``2**24`` of ``geometric``, ``logseries`` and ``zipf``, and the ones of
+  ``poisson`` and ``binomial`` with a standard deviation above about
+  ``2**22``.
+* ``zipf`` follows the Zipf distribution truncated to the largest ``int64``.
+  NumPy's rejection test rounds to accepting every sample above ``2**53``,
+  which makes its tail too heavy for ``a`` close to 1.
+* :meth:`RandomState.randint` and :func:`randint` return ``int64`` values by
+  default, as do :meth:`RandomState.random_integers`,
+  :func:`random_integers` and :meth:`RandomState.tomaxint`, whose values
+  reach the largest ``int64``. NumPy returns a C ``long`` instead, which is
+  ``int32`` on Windows.
+* ``hypergeometric`` requires ``ngood`` and ``nbad`` below ``10**9``, also in
+  :class:`RandomState`.
+* :meth:`Generator.shuffle` and :func:`shuffle` only accept a
+  :class:`dpnp.ndarray`, as in CuPy.
+
+Random Generator
+----------------
+
+.. autosummary::
+   :toctree: generated/
+   :nosignatures:
+
+   default_rng
+   Generator
+
+Simple random data
+^^^^^^^^^^^^^^^^^^
+
+.. autosummary::
+   :toctree: generated/
+   :nosignatures:
+
+   Generator.integers
+   Generator.random
+   Generator.choice
+   Generator.bytes
+
+Permutations
+^^^^^^^^^^^^
+
+.. autosummary::
+   :toctree: generated/
+   :nosignatures:
+
+   Generator.shuffle
+   Generator.permutation
+   Generator.permuted
+
+Distributions
+^^^^^^^^^^^^^
+
+.. autosummary::
+   :toctree: generated/
+   :nosignatures:
+
+   Generator.beta
+   Generator.binomial
+   Generator.chisquare
+   Generator.dirichlet
+   Generator.exponential
+   Generator.f
+   Generator.gamma
+   Generator.geometric
+   Generator.gumbel
+   Generator.hypergeometric
+   Generator.laplace
+   Generator.logistic
+   Generator.lognormal
+   Generator.logseries
+   Generator.multinomial
+   Generator.multivariate_hypergeometric
+   Generator.multivariate_normal
+   Generator.negative_binomial
+   Generator.noncentral_chisquare
+   Generator.noncentral_f
+   Generator.normal
+   Generator.pareto
+   Generator.poisson
+   Generator.power
+   Generator.rayleigh
+   Generator.standard_cauchy
+   Generator.standard_exponential
+   Generator.standard_gamma
+   Generator.standard_normal
+   Generator.standard_t
+   Generator.triangular
+   Generator.uniform
+   Generator.vonmises
+   Generator.wald
+   Generator.weibull
+   Generator.zipf
+
+Bit generators
+--------------
+
+.. autosummary::
+   :toctree: generated/
+   :nosignatures:
+
+   BitGenerator
+   Philox4x3210
+   PCG64DXSM
+   MRG32k3a
+   SeedSequence
+
+Legacy random generation
+------------------------
+
 .. hint:: `NumPy API Reference: Legacy random generation <https://numpy.org/doc/stable/reference/random/legacy.html>`_
 
+:class:`RandomState` keeps the parameter checks of NumPy's legacy one, also
+for parameters in :class:`dpnp.ndarray`, as CuPy's does, so that a call with
+such parameters waits for the device to check them. The module-level
+functions use one :class:`RandomState` per SYCL queue: the
+queue given by the ``device`` or ``sycl_queue`` keyword, or else the queue of
+a :class:`dpnp.ndarray` parameter, or else the default queue.
+
+.. autosummary::
+   :toctree: generated/
+   :nosignatures:
+
+   RandomState
+   seed
+   get_state
+   set_state
+   get_bit_generator
+   set_bit_generator
+
 Simple random data
-------------------
+^^^^^^^^^^^^^^^^^^
 
 .. autosummary::
    :toctree: generated/
@@ -23,9 +198,8 @@ Simple random data
    choice
    bytes
 
-
 Permutations
-------------
+^^^^^^^^^^^^
 
 .. autosummary::
    :toctree: generated/
@@ -34,9 +208,8 @@ Permutations
    shuffle
    permutation
 
-
 Distributions
--------------
+^^^^^^^^^^^^^
 
 .. autosummary::
    :toctree: generated/
@@ -77,16 +250,3 @@ Distributions
    wald
    weibull
    zipf
-
-
-Random generator
-----------------
-
-.. autosummary::
-   :toctree: generated/
-   :nosignatures:
-
-   RandomState
-   seed
-   get_random_state
-   set_random_state

dpnp/AGENTS.md

@@ -6,7 +6,8 @@
 - `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`
-- 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`).
+- `random/`: `_bit_generator.py` (`BitGenerator`, `Philox4x3210`, `PCG64DXSM`, `MRG32k3a`; `SeedSequence` is NumPy's), `_generator.py` (`Generator`, which implements every distribution), `dpnp_random_state.py` (`RandomState`, the same samplers with NumPy's legacy checks) and `dpnp_iface_random.py` (module functions, one `RandomState` per queue). 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`). New kernels go into `dpnp/libtensor` (see `dpnp/libtensor/AGENTS.md`).
 
 ## Adding a Public Function
 1. Implement it in the matching `dpnp_iface_*.py` or subpackage module. Iface modules have no `__all__`.
@@ -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`, `_float_sample` or `_discrete_sample` into `BitGenerator._fill`. `RandomState` and the module function delegate to it. Parameters on the device are checked only by `RandomState` (`_check_device`), as the check waits for the device. Temporaries take the output's `usm_type`.
 - 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,9 +229,9 @@ 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)
 
 add_subdirectory(dpnp_algo)
 add_subdirectory(dpnp_utils)
-add_subdirectory(random)

dpnp/backend/AGENTS.md

@@ -3,5 +3,5 @@
 Legacy, don't extend: new kernels and bindings go into `dpnp/libtensor` (see `dpnp/libtensor/AGENTS.md`, which also holds the C++ rules and build options).
 
 ## Layout
-- `src/`, `include/` and `kernels/dpnp_krnl_*.cpp` build `dpnp_backend_c` (still used by `dpnp.partition` and `dpnp.random`); `tests/` (GTest, `-DDPNP_BACKEND_TESTS=ON`) and `examples/` are not built in CI
+- `src/`, `include/` and `kernels/dpnp_krnl_*.cpp` build `dpnp_backend_c` (still used by `dpnp.partition`)
 - `cmake/Modules/`: vendored `IntelSYCL`, `MKL`, `TBB` and `oneDPL` config modules, used by the whole build

dpnp/backend/CMakeLists.txt

@@ -27,17 +27,7 @@
 # THE POSSIBILITY OF SUCH DAMAGE.
 # *****************************************************************************
 
-set(DPNP_SRC
-    kernels/dpnp_krnl_arraycreation.cpp
-    kernels/dpnp_krnl_common.cpp
-    kernels/dpnp_krnl_random.cpp
-    kernels/dpnp_krnl_sorting.cpp
-    src/dpnp_iface_fptr.cpp
-    src/memory_sycl.cpp
-    src/queue_sycl.cpp
-    src/verbose.cpp
-    src/dpnp_random_state.cpp
-)
+set(DPNP_SRC kernels/dpnp_krnl_sorting.cpp src/dpnp_iface_fptr.cpp)
 
 set(_trgt "dpnp_backend_c")
 
@@ -69,12 +59,6 @@ if(DPNP_GENERATE_COVERAGE)
     target_link_options(${_trgt} PRIVATE -fprofile-instr-generate -fcoverage-mapping)
 endif()
 
-# Link only the MKL SYCL domains actually used by the backend sources
-# (BLAS: mkl_blas::scal; RNG: engines/distributions; VM: vector math ops)
-target_link_libraries(
-    ${_trgt}
-    PUBLIC MKL::MKL_SYCL::BLAS MKL::MKL_SYCL::RNG MKL::MKL_SYCL::VM
-)
 target_link_libraries(${_trgt} PUBLIC oneDPL)
 
 if(UNIX)
@@ -103,10 +87,6 @@ if(DPNP_WITH_REDIST)
     set_target_properties(${_trgt} PROPERTIES INSTALL_RPATH "$ORIGIN/../../../")
 endif()
 
-if(DPNP_BACKEND_TESTS)
-    add_subdirectory(tests)
-endif()
-
 set(_lib_destination dpnp)
 set(_include_destination dpnp/backend/include)
 install(