Add PCG64DXSM and MRG32k3a bit generators to dpnp.random

CHANGELOG.md

@@ -14,6 +14,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 * Added the continuous distributions of `numpy.random.Generator`, such as `gamma`, `beta` and `vonmises`, to `dpnp.random.Generator` [#3158](https://github.com/IntelPython/dpnp/pull/3158)
 * Added `integers` and the discrete distributions of `numpy.random.Generator`, such as `poisson`, `binomial` and `hypergeometric`, to `dpnp.random.Generator`, along with `noncentral_chisquare` and `noncentral_f` [#3159](https://github.com/IntelPython/dpnp/pull/3159)
 * Added the multivariate distributions, `permutation`, `shuffle`, `permuted`, `choice` and `bytes` to `dpnp.random.Generator` [#3160](https://github.com/IntelPython/dpnp/pull/3160)
+* Added the `dpnp.random.PCG64DXSM` bit generator, which matches `numpy.random.PCG64DXSM` for raw values, `random`, `uniform` and `bytes` [#3161](https://github.com/IntelPython/dpnp/pull/3161)
+* Added the `dpnp.random.MRG32k3a` bit generator [#3161](https://github.com/IntelPython/dpnp/pull/3161)
 
 ### Changed
 

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

@@ -5,8 +5,10 @@ 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. It provides the
+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 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.
@@ -31,8 +33,10 @@ 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.
+  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.
@@ -144,6 +148,8 @@ Bit generators
 
    BitGenerator
    Philox4x3210
+   PCG64DXSM
+   MRG32k3a
    SeedSequence
 
 Legacy random generation

dpnp/AGENTS.md

@@ -6,7 +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`.
+- `random/`: `_bit_generator.py` (`BitGenerator`, `Philox4x3210`, `PCG64DXSM`, `MRG32k3a`; `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

dpnp/libtensor/extensions/rng/common.hpp

@@ -128,7 +128,7 @@ sycl::event submit_sample(sycl::queue &exec_q,
                           const std::size_t nelems,
                           const std::vector<sycl::event> &depends)
 {
-    const std::size_t n_wi = rng_ns::num_work_items<EngineTraitsT>(nelems);
+    const std::size_t n_wi = rng_ns::num_work_items<EngineTraitsT>(nelems, L);
 
     return exec_q.submit([&](sycl::handler &cgh) {
         cgh.depends_on(depends);
@@ -182,7 +182,12 @@ sycl::event sample_impl(sycl::queue &exec_q,
 
     T *res = reinterpret_cast<T *>(result);
 
-    if constexpr (rng_ns::can_draw_blocked<EngineTraitsT, SamplerT>) {
+    if constexpr (rng_ns::can_draw_sequential<EngineTraitsT, SamplerT>) {
+        return detail::submit_sample<EngineTraitsT, SamplerT,
+                                     rng_ns::Layout::sequential>(
+            exec_q, state, sampler, res, nelems, depends);
+    }
+    else if constexpr (rng_ns::can_draw_blocked<EngineTraitsT, SamplerT>) {
         // the blocked layout gives the same values, faster on GPUs but
         // slower on CPUs
         if (exec_q.get_device().is_gpu()) {
@@ -242,26 +247,48 @@ void init_sample_dispatch_table(sample_dispatch_table_t &table)
 {
     using ext::common::init_dispatch_vector;
 
-    static_assert(rng_ns::num_engines == 1,
+    static_assert(rng_ns::num_engines == 3,
                   "Every engine must have a row in the dispatch table");
 
     constexpr int philox = static_cast<int>(rng_ns::engine_id::philox4x32x10);
     init_dispatch_vector<sample_fn_ptr_t,
                          SampleFactoryFor<rng_ns::Philox4x32x10, SamplerT,
                                           SupportedT>::template type>(
         table[philox]);
+
+    constexpr int pcg = static_cast<int>(rng_ns::engine_id::pcg64dxsm);
+    init_dispatch_vector<sample_fn_ptr_t,
+                         SampleFactoryFor<rng_ns::PCG64DXSM, SamplerT,
+                                          SupportedT>::template type>(
+        table[pcg]);
+
+    constexpr int mrg = static_cast<int>(rng_ns::engine_id::mrg32k3a);
+    init_dispatch_vector<sample_fn_ptr_t,
+                         SampleFactoryFor<rng_ns::MRG32k3a, SamplerT,
+                                          SupportedT>::template type>(
+        table[mrg]);
 }
 
 template <template <typename> class SamplerT>
 void init_discrete_dispatch_table(discrete_dispatch_table_t &table)
 {
-    static_assert(rng_ns::num_engines == 1,
+    static_assert(rng_ns::num_engines == 3,
                   "Every engine must have a row in the dispatch table");
 
     using Philox = rng_ns::Philox4x32x10;
     constexpr int philox = static_cast<int>(rng_ns::engine_id::philox4x32x10);
     table[philox][0] = sample_impl<Philox, SamplerT<float>, float>;
     table[philox][1] = sample_impl<Philox, SamplerT<double>, double>;
+
+    using PCG = rng_ns::PCG64DXSM;
+    constexpr int pcg = static_cast<int>(rng_ns::engine_id::pcg64dxsm);
+    table[pcg][0] = sample_impl<PCG, SamplerT<float>, float>;
+    table[pcg][1] = sample_impl<PCG, SamplerT<double>, double>;
+
+    using MRG = rng_ns::MRG32k3a;
+    constexpr int mrg = static_cast<int>(rng_ns::engine_id::mrg32k3a);
+    table[mrg][0] = sample_impl<MRG, SamplerT<float>, float>;
+    table[mrg][1] = sample_impl<MRG, SamplerT<double>, double>;
 }
 
 // Looks up the sampler of `result`'s type, and throws if there is none