From c5182f4050f83d79884182169a0d21769881b379 Mon Sep 17 00:00:00 2001
From: Nikita Grigorian <nikita.grigorian@intel.com>
Date: Wed, 4 Nov 2026 11:46:56 -0800
Subject: [PATCH] Add the PCG64DXSM and MRG32k3a bit generators to dpnp.random

PCG64DXSM is a custom engine bit-identical with NumPy's for raw values,
random, uniform and bytes, drawn in a sequential layout for the samplers
taking a fixed number of draws. MRG32k3a uses oneMKL's device engine,
with work-item w on substream w.

sycl::mul_hi is replaced with a limb-based umul_hi, as DPC++ 2026.2
miscompiles it on Intel GPUs.
---
 CHANGELOG.md                                  |   2 +
 doc/known_words.txt                           |   6 +
 doc/reference/random.rst                      |  21 +-
 dpnp/AGENTS.md                                |   2 +-
 dpnp/libtensor/extensions/rng/common.hpp      |  35 +-
 .../extensions/rng/kernels/engines.hpp        | 236 +++++++++-
 .../extensions/rng/kernels/sample.hpp         |  45 +-
 .../extensions/rng/kernels/samplers.hpp       |  11 +-
 dpnp/libtensor/extensions/rng/rng_py.cpp      |   8 +-
 dpnp/random/__init__.py                       |  10 +-
 dpnp/random/_bit_generator.py                 | 407 ++++++++++++++++++
 dpnp/random/_generator.py                     |   8 +-
 dpnp/tests/AGENTS.md                          |   2 +-
 dpnp/tests/test_random_generator.py           | 346 ++++++++++++++-
 dpnp/tests/test_sycl_queue.py                 |   5 +-
 dpnp/tests/test_usm_type.py                   |   5 +-
 .../cupy/random_tests/test_bit_generator.py   |   1 -
 17 files changed, 1118 insertions(+), 32 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index bcaf298232c..6bc0ed02c78 100644
--- a/CHANGELOG.md
+++ b/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
 
diff --git a/doc/known_words.txt b/doc/known_words.txt
index 4d81cd5b853..d537c6ef46b 100644
--- a/doc/known_words.txt
+++ b/doc/known_words.txt
@@ -43,9 +43,11 @@ Hanning
 histogrammed
 Hypergeometric
 kwargs
+Kelton
 iaxis
 iinfo
 Infs
+L'Ecuyer
 intp
 ints
 iterable
@@ -89,17 +91,21 @@ prepending
 reentrant
 representable
 resampling
+RngStreams
 runtimes
 scikit
 se
 SeedSequence
 signbit
 signum
+Simard
 sinc
 strided
 subarray
 subarrays
 subclasses
+substream
+substreams
 subtype
 SyclDevice
 SyclQueue
diff --git a/doc/reference/random.rst b/doc/reference/random.rst
index 93da08ce480..a937481ded1 100644
--- a/doc/reference/random.rst
+++ b/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.
@@ -28,7 +30,10 @@ A generator draws all its samples on one SYCL queue, chosen with the
 Differences from NumPy
 ----------------------
 
-* The streams differ from NumPy's. For a given seed, output size and data
+* 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`. For a given seed, output size and data
   type, the raw bits and the uniform floats are the same on every device.
   The other distributions use transcendental functions, which may round
   differently on another device: their samples agree to a few ULP, or
@@ -37,7 +42,8 @@ Differences from NumPy
   precision, as noted below.
 * 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.
+  ``2 * n``, but for the methods :class:`PCG64DXSM` draws in the order of
+  its stream. 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
@@ -66,8 +72,9 @@ Differences from NumPy
 * ``hypergeometric`` requires ``ngood`` and ``nbad`` below ``10**9``.
 * :meth:`Generator.shuffle` only accepts a :class:`dpnp.ndarray`, as in
   CuPy.
-* Generators and bit generators cannot be pickled or copied, and the bit
-  generators have no ``advance`` or ``jumped`` method.
+* Generators and bit generators cannot be pickled or copied, and
+  :class:`Philox4x3210` and :class:`MRG32k3a` have no ``advance`` or
+  ``jumped`` method.
 
 Random Generator
 ----------------
@@ -168,6 +175,8 @@ Bit generators
    BitGenerator.spawn
    BitGenerator.state
    Philox4x3210
+   PCG64DXSM
+   MRG32k3a
    SeedSequence
 
 Legacy random generation
diff --git a/dpnp/AGENTS.md b/dpnp/AGENTS.md
index 6009def4f08..3f545cb7ef4 100644
--- a/dpnp/AGENTS.md
+++ b/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
diff --git a/dpnp/libtensor/extensions/rng/common.hpp b/dpnp/libtensor/extensions/rng/common.hpp
index abdd80042eb..4c439680748 100644
--- a/dpnp/libtensor/extensions/rng/common.hpp
+++ b/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,7 +247,7 @@ 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);
@@ -250,18 +255,40 @@ void init_sample_dispatch_table(sample_dispatch_table_t &table)
                          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
diff --git a/dpnp/libtensor/extensions/rng/kernels/engines.hpp b/dpnp/libtensor/extensions/rng/kernels/engines.hpp
index c9af142bd7d..bc7791e9226 100644
--- a/dpnp/libtensor/extensions/rng/kernels/engines.hpp
+++ b/dpnp/libtensor/extensions/rng/kernels/engines.hpp
@@ -36,12 +36,14 @@
 
 #include <oneapi/mkl/rng/device.hpp>
 
+#include "umul_hi.hpp"
+
 namespace dpnp::kernels::rng
 {
 namespace mkl_rng_dev = oneapi::mkl::rng::device;
 
 // The maximum number of 64-bit seed words any supported engine takes.
-inline constexpr std::size_t max_seed_words = 1;
+inline constexpr std::size_t max_seed_words = 10;
 
 // Host-side state of a bit generator, passed by value into kernels.
 //
@@ -51,7 +53,8 @@ inline constexpr std::size_t max_seed_words = 1;
 // other work-item, of this or any later launch, no matter how many draws a
 // rejection sampler takes. For Philox, work-item `w` of a launch owns the
 // counter block `counter + w * 2**64`, in units of 128-bit counter
-// increments.
+// increments. PCG64DXSM and MRG32k3a keep their whole state in the seed
+// words.
 struct EngineState
 {
     std::uint64_t seed[max_seed_words];
@@ -64,8 +67,10 @@ struct EngineState
 enum class engine_id : int
 {
     philox4x32x10 = 0,
+    pcg64dxsm = 1,
+    mrg32k3a = 2,
 };
-inline constexpr int num_engines = 1;
+inline constexpr int num_engines = 3;
 
 // The conversions of raw bits to uniform floating-point values, for a bit
 // generator `Derived` providing next_uint32() and next_uint64(). They are
@@ -170,6 +175,10 @@ struct Philox4x32x10
     using engine_t = mkl_rng_dev::philox4x32x10<1>;
     using bitgen_t = BitGen<engine_t>;
 
+    // Whether the fixed-consumption samplers draw the stream in order,
+    // see SampleFunctor
+    static constexpr bool ordered = false;
+
     // The elements a work-item produces in the strided and blocked layouts
     static constexpr std::size_t items_per_work_item = 8;
 
@@ -199,4 +208,225 @@ struct Philox4x32x10
     }
 };
 
+// An unsigned 128-bit integer, with the arithmetic modulo 2**128 the LCG of
+// PCG64 needs
+struct UInt128
+{
+    std::uint64_t lo;
+    std::uint64_t hi;
+};
+
+inline UInt128 operator+(const UInt128 &a, const UInt128 &b)
+{
+    const std::uint64_t lo = a.lo + b.lo;
+    return {lo, a.hi + b.hi + (lo < a.lo ? 1 : 0)};
+}
+
+inline UInt128 operator*(const UInt128 &a, const UInt128 &b)
+{
+    return {a.lo * b.lo, umul_hi(a.lo, b.lo) + a.lo * b.hi + a.hi * b.lo};
+}
+
+// The affine map x -> mult * x + inc modulo 2**128, as a step of an LCG
+struct Affine128
+{
+    UInt128 mult;
+    UInt128 inc;
+
+    UInt128 operator()(const UInt128 &x) const { return mult * x + inc; }
+};
+
+// The map `m` applied `k` times, in O(log k) steps (F. Brown, "Random
+// number generation with arbitrary stride", 1994)
+inline Affine128 affine_power(Affine128 m, std::uint64_t k)
+{
+    Affine128 acc{{1, 0}, {0, 0}};
+    while (k != 0) {
+        if (k & 1) {
+            acc = {acc.mult * m.mult, acc.inc * m.mult + m.inc};
+        }
+        m = {m.mult * m.mult, (m.mult + UInt128{1, 0}) * m.inc};
+        k >>= 1;
+    }
+    return acc;
+}
+
+// NumPy's PCG64DXSM: a 128-bit LCG with the cheap 64-bit multiplier and an
+// odd increment, whose output is the DXSM permutation of the state before
+// the step. 32-bit draws take the low half of an output, and keep its high
+// half for the next one, as NumPy does.
+class Pcg64Dxsm : public UniformBits<Pcg64Dxsm>
+{
+private:
+    UInt128 state_;
+    UInt128 inc_;
+    bool has_uint32_ = false;
+    std::uint32_t uinteger_ = 0;
+
+public:
+    static constexpr std::uint64_t cheap_multiplier = 0xda942042e4dd58b5ULL;
+
+    Pcg64Dxsm(const UInt128 &state, const UInt128 &inc)
+        : state_(state), inc_(inc)
+    {
+    }
+
+    // Use `value` for the next 32-bit draw
+    void set_buffered(std::uint32_t value)
+    {
+        has_uint32_ = true;
+        uinteger_ = value;
+    }
+
+    std::uint64_t next_uint64()
+    {
+        std::uint64_t hi = state_.hi;
+        const std::uint64_t lo = state_.lo | 1;
+        state_ = state_ * UInt128{cheap_multiplier, 0} + inc_;
+
+        hi ^= hi >> 32;
+        hi *= cheap_multiplier;
+        hi ^= hi >> 48;
+        hi *= lo;
+        return hi;
+    }
+
+    std::uint32_t next_uint32()
+    {
+        if (has_uint32_) {
+            has_uint32_ = false;
+            return uinteger_;
+        }
+        const std::uint64_t next = next_uint64();
+        set_buffered(static_cast<std::uint32_t>(next >> 32));
+        return static_cast<std::uint32_t>(next);
+    }
+};
+
+struct PCG64DXSM
+{
+    using bitgen_t = Pcg64Dxsm;
+
+    // The samplers that take a fixed number of 32-bit halves of raw outputs
+    // per element draw the stream in order, as NumPy does
+    static constexpr bool ordered = true;
+
+    // The elements a work-item produces in the strided layout
+    static constexpr std::size_t items_per_work_item = 8;
+
+    // No make_block
+    static constexpr unsigned max_block_words = 0;
+
+    // The seed words are the state, the increment, the map of the stride of
+    // NumPy's jumped() (the step applied (phi - 1) * 2**128 times), each as
+    // {lo, hi}, then NumPy's has_uint32 and uinteger.
+    static UInt128 word128(const EngineState &state, std::size_t i)
+    {
+        return {state.seed[i], state.seed[i + 1]};
+    }
+
+    // The stream jumped `subsequence` times, as by NumPy's jumped(), for the
+    // samplers whose draws per element vary
+    static bitgen_t make(const EngineState &state, std::uint64_t subsequence)
+    {
+        const Affine128 stride{word128(state, 4), word128(state, 6)};
+        return bitgen_t(affine_power(stride, subsequence)(word128(state, 0)),
+                        word128(state, 2));
+    }
+
+    // The stream at element `first` of a sampler taking `halves` 32-bit
+    // halves of raw outputs per element
+    static bitgen_t
+        make_at(const EngineState &state, std::uint64_t first, unsigned halves)
+    {
+        const UInt128 inc = word128(state, 2);
+        const Affine128 step{{Pcg64Dxsm::cheap_multiplier, 0}, inc};
+        const UInt128 start = word128(state, 0);
+
+        if (halves == 2) {
+            // 64-bit draws, which leave the buffered half alone
+            return bitgen_t(affine_power(step, first)(start), inc);
+        }
+
+        // 32-bit draws, the first of which is the buffered half, if any
+        const bool has_uint32 = state.seed[8] != 0;
+        if (has_uint32 && first == 0) {
+            bitgen_t bg(start, inc);
+            bg.set_buffered(static_cast<std::uint32_t>(state.seed[9]));
+            return bg;
+        }
+        const std::uint64_t half = has_uint32 ? first - 1 : first;
+        bitgen_t bg(affine_power(step, half / 2)(start), inc);
+        if (half % 2 != 0) {
+            // start with the high half of an output
+            bg.next_uint32();
+        }
+        return bg;
+    }
+};
+
+// Adapts oneMKL's mrg32k3a, whose raw values r are in [1, m1] with
+// m1 = 2**32 - 209, to uniform 32-bit words. v = r - 1 is spread over the
+// 32-bit range as v + floor(209 * v / 2**32), which increases strictly, so
+// every word has probability 1 / m1, except 209 words that never occur.
+template <typename EngineT>
+class Mrg32k3aBitGen : public UniformBits<Mrg32k3aBitGen<EngineT>>
+{
+private:
+    // 2**32 - m1
+    static constexpr std::uint64_t gap = 209;
+
+    EngineT engine_;
+    mkl_rng_dev::bits<std::uint32_t> bits_{};
+
+public:
+    explicit Mrg32k3aBitGen(const EngineT &engine) : engine_(engine) {}
+
+    std::uint32_t next_uint32()
+    {
+        const std::uint64_t v = mkl_rng_dev::generate(bits_, engine_) - 1;
+        return static_cast<std::uint32_t>(v + ((v * gap) >> 32));
+    }
+
+    std::uint64_t next_uint64()
+    {
+        const std::uint64_t hi = next_uint32();
+        const std::uint64_t lo = next_uint32();
+        return (hi << 32) | lo;
+    }
+};
+
+// L'Ecuyer's MRG32k3a combined multiple recursive generator, as implemented
+// by the oneMKL device RNG API
+struct MRG32k3a
+{
+    using engine_t = mkl_rng_dev::mrg32k3a<1>;
+    using bitgen_t = Mrg32k3aBitGen<engine_t>;
+
+    static constexpr bool ordered = false;
+
+    // The elements a work-item produces in the strided layout. Moving the
+    // engine to the substream of a work-item takes a 3x3 matrix product
+    // modulo each of m1 and m2 per base-8 digit of the work-item's index,
+    // which 8 elements per work-item made 2x (GPU) to 10x (CPU) slower than
+    // Philox. 128 is on par with Philox on both.
+    static constexpr std::size_t items_per_work_item = 128;
+
+    // No make_block: the engine computes its values one after another
+    static constexpr unsigned max_block_words = 0;
+
+    // The seed words are the state at the start of the launch, the three
+    // words of the first component, then those of the second. Work-item
+    // `w` draws from substream `w`, which starts 2**76 * w values later, as
+    // the substreams of L'Ecuyer et al.'s RngStreams.
+    static bitgen_t make(const EngineState &state, std::uint64_t subsequence)
+    {
+        const auto word = [&state](std::size_t i) {
+            return static_cast<std::uint32_t>(state.seed[i]);
+        };
+        return bitgen_t(
+            engine_t({word(0), word(1), word(2), word(3), word(4), word(5)},
+                     {std::uint64_t{0}, subsequence << 12, subsequence >> 52}));
+    }
+};
 } // namespace dpnp::kernels::rng
diff --git a/dpnp/libtensor/extensions/rng/kernels/sample.hpp b/dpnp/libtensor/extensions/rng/kernels/sample.hpp
index 1f99690cc34..58279b0294e 100644
--- a/dpnp/libtensor/extensions/rng/kernels/sample.hpp
+++ b/dpnp/libtensor/extensions/rng/kernels/sample.hpp
@@ -36,6 +36,11 @@
 
 namespace dpnp::kernels::rng
 {
+// The number of output elements each work-item produces in the sequential
+// layout, where it only sets how much of the stream every work-item skips
+// to, not the values.
+inline constexpr std::size_t sequential_items_per_work_item = 64;
+
 // How the elements of a launch map to engine streams:
 // - strided: work-item `w` builds its own engine on the counter block `w`
 //   and draws elements w, w + G, w + 2G, ..., where G is the number of
@@ -45,7 +50,12 @@ namespace dpnp::kernels::rng
 //   of the stream definition: changing it changes every generated value;
 // - blocked: the strided layout, with all the words of a work-item
 //   generated together, which is faster on GPUs. It needs a sampler taking
-//   a fixed number of words per element, and gives the strided values.
+//   a fixed number of words per element, and gives the strided values;
+// - sequential: for an ordered engine and a sampler taking a fixed number
+//   of draws per element, the elements are the sequential stream of the
+//   engine, and work-item `w` draws elements w * S to (w + 1) * S - 1,
+//   where S is sequential_items_per_work_item, from the engine advanced
+//   past the earlier ones.
 // The words each element draws depend only on the engine state and the
 // number of elements, never on the device, the layout chosen for it or how
 // the work-items are scheduled. Raw bits and uniform floats are therefore
@@ -56,8 +66,13 @@ enum class Layout
 {
     strided,
     blocked,
+    sequential,
 };
 
+template <typename EngineTraitsT, typename SamplerT>
+inline constexpr bool can_draw_sequential =
+    EngineTraitsT::ordered && SamplerT::halves_per_item != 0;
+
 template <typename EngineTraitsT, typename SamplerT>
 inline constexpr bool can_draw_blocked =
     SamplerT::words_per_item != 0 &&
@@ -66,9 +81,11 @@ inline constexpr bool can_draw_blocked =
 
 // The number of work-items a launch producing `n` elements uses
 template <typename EngineTraitsT>
-inline std::size_t num_work_items(std::size_t n)
+inline std::size_t num_work_items(std::size_t n, Layout layout)
 {
-    constexpr std::size_t m = EngineTraitsT::items_per_work_item;
+    const std::size_t m = (layout == Layout::sequential)
+                              ? sequential_items_per_work_item
+                              : EngineTraitsT::items_per_work_item;
     return (n + m - 1) / m;
 }
 
@@ -79,6 +96,8 @@ class SampleFunctor
 private:
     using T = typename SamplerT::result_type;
 
+    static_assert(L != Layout::sequential ||
+                  can_draw_sequential<EngineTraitsT, SamplerT>);
     static_assert(L != Layout::blocked ||
                   can_draw_blocked<EngineTraitsT, SamplerT>);
 
@@ -106,7 +125,25 @@ class SampleFunctor
         const std::size_t w = id[0];
         SamplerT sampler = sampler_;
 
-        if constexpr (L == Layout::blocked) {
+        if constexpr (L == Layout::sequential) {
+            constexpr std::size_t m = sequential_items_per_work_item;
+            const std::size_t first = w * m;
+
+            auto bg = EngineTraitsT::make_at(state_, first,
+                                             SamplerT::halves_per_item);
+            if (first + m <= n_) {
+#pragma unroll
+                for (std::size_t k = 0; k < m; ++k) {
+                    out_[first + k] = sampler(bg, first + k);
+                }
+            }
+            else {
+                for (std::size_t i = first; i < n_; ++i) {
+                    out_[i] = sampler(bg, i);
+                }
+            }
+        }
+        else if constexpr (L == Layout::blocked) {
             constexpr std::size_t items = EngineTraitsT::items_per_work_item;
             constexpr unsigned words = SamplerT::words_per_item * items;
 
diff --git a/dpnp/libtensor/extensions/rng/kernels/samplers.hpp b/dpnp/libtensor/extensions/rng/kernels/samplers.hpp
index 6d651d324a4..93ac6e6aa7c 100644
--- a/dpnp/libtensor/extensions/rng/kernels/samplers.hpp
+++ b/dpnp/libtensor/extensions/rng/kernels/samplers.hpp
@@ -46,6 +46,10 @@ namespace dpnp::kernels::rng
 // `result_type` is the type of the samples, and `param_types` lists the
 // type of each parameter, in the order the constructor takes them.
 //
+// `halves_per_item` is the number of 32-bit halves of raw 64-bit outputs
+// the sampler takes for every element, if that number is fixed, and 0
+// otherwise.
+//
 // `words_per_item` is the number of 32-bit words the sampler takes per
 // element, on average over the items_per_work_item consecutive elements of one
 // work-item, if that number is fixed whatever the parameters and the
@@ -57,6 +61,7 @@ struct SamplerBase
     using result_type = ResT;
     using param_types = std::tuple<ParamTs...>;
 
+    static constexpr unsigned halves_per_item = 0;
     static constexpr unsigned words_per_item = 0;
 };
 
@@ -69,7 +74,8 @@ template <typename T>
 class RawSampler : public SamplerBase<T>
 {
 public:
-    static constexpr unsigned words_per_item = sizeof(T) / 4;
+    static constexpr unsigned halves_per_item = sizeof(T) / 4;
+    static constexpr unsigned words_per_item = halves_per_item;
 
     template <typename BitGenT>
     T operator()(BitGenT &bg, const std::size_t)
@@ -92,7 +98,8 @@ class UniformSampler : public SamplerBase<T, T, T>
     Param<T> high_;
 
 public:
-    static constexpr unsigned words_per_item = sizeof(T) / 4;
+    static constexpr unsigned halves_per_item = sizeof(T) / 4;
+    static constexpr unsigned words_per_item = halves_per_item;
 
     UniformSampler(const Param<T> &low, const Param<T> &high)
         : low_(low), high_(high)
diff --git a/dpnp/libtensor/extensions/rng/rng_py.cpp b/dpnp/libtensor/extensions/rng/rng_py.cpp
index 512243045ec..0e85998bc48 100644
--- a/dpnp/libtensor/extensions/rng/rng_py.cpp
+++ b/dpnp/libtensor/extensions/rng/rng_py.cpp
@@ -45,13 +45,19 @@ PYBIND11_MODULE(_rng_impl, m, py::mod_gil_not_used())
     using dpnp::kernels::rng::engine_id;
 
     py::enum_<engine_id>(m, "engine_id")
-        .value("philox4x32x10", engine_id::philox4x32x10);
+        .value("philox4x32x10", engine_id::philox4x32x10)
+        .value("pcg64dxsm", engine_id::pcg64dxsm)
+        .value("mrg32k3a", engine_id::mrg32k3a);
 
     // The elements a work-item produces in the strided layout, by engine
     namespace kr = dpnp::kernels::rng;
     py::dict items_per_work_item;
     items_per_work_item[py::cast(engine_id::philox4x32x10)] =
         kr::Philox4x32x10::items_per_work_item;
+    items_per_work_item[py::cast(engine_id::pcg64dxsm)] =
+        kr::PCG64DXSM::items_per_work_item;
+    items_per_work_item[py::cast(engine_id::mrg32k3a)] =
+        kr::MRG32k3a::items_per_work_item;
     m.attr("items_per_work_item") = items_per_work_item;
 
     rng_ns::init_uniform(m);
diff --git a/dpnp/random/__init__.py b/dpnp/random/__init__.py
index 4e7a50ad44d..dbf86094885 100644
--- a/dpnp/random/__init__.py
+++ b/dpnp/random/__init__.py
@@ -26,7 +26,13 @@
 # THE POSSIBILITY OF SUCH DAMAGE.
 # *****************************************************************************
 
-from ._bit_generator import BitGenerator, Philox4x3210, SeedSequence
+from ._bit_generator import (
+    PCG64DXSM,
+    BitGenerator,
+    MRG32k3a,
+    Philox4x3210,
+    SeedSequence,
+)
 from ._generator import Generator, default_rng
 from .dpnp_iface_random import (
     beta,
@@ -83,6 +89,8 @@
 __all__ = [
     "BitGenerator",
     "Generator",
+    "MRG32k3a",
+    "PCG64DXSM",
     "Philox4x3210",
     "SeedSequence",
     "default_rng",
diff --git a/dpnp/random/_bit_generator.py b/dpnp/random/_bit_generator.py
index 479b9b84630..95b19fb6ac5 100644
--- a/dpnp/random/_bit_generator.py
+++ b/dpnp/random/_bit_generator.py
@@ -40,6 +40,7 @@
 # pylint: disable=no-name-in-module
 # pylint: disable=protected-access
 
+import functools
 import math
 import operator
 import threading
@@ -52,6 +53,8 @@
 
 __all__ = [
     "BitGenerator",
+    "MRG32k3a",
+    "PCG64DXSM",
     "Philox4x3210",
     "SeedSequence",
 ]
@@ -447,3 +450,407 @@ def _launch(self, kernel, dtype, size):
             counter = self._counter
             self._counter = (counter + (n_blocks << 64)) & _MASK128
         return [key], counter & _MASK64, counter >> 64
+
+
+# NumPy's PCG64DXSM: the multiplier of its LCG, and the step of jumped(),
+# (phi - 1) * 2**128 made odd
+_PCG_MULT = 0xDA942042E4DD58B5
+_PCG_JUMP = 0x9E3779B97F4A7C15F39CC0605CEDC835
+
+# The kernels that take a fixed number of 32-bit halves of raw outputs for
+# every value of a data type, and so draw the stream of PCG64DXSM in order,
+# as NumPy does. It must match halves_per_item of the samplers.
+_PCG_SEQUENTIAL = {
+    ri._random_raw: {
+        numpy.dtype(numpy.uint64): 2,
+        numpy.dtype(numpy.uint32): 1,
+    },
+    ri._uniform: {numpy.dtype(numpy.float64): 2, numpy.dtype(numpy.float32): 1},
+}
+
+
+def _pcg_power(inc, k):
+    """The step of the LCG applied `k` times, as a map (mult, inc)."""
+
+    acc_mult, acc_inc = 1, 0
+    mult = _PCG_MULT
+    while k:
+        if k & 1:
+            acc_mult = (acc_mult * mult) & _MASK128
+            acc_inc = (acc_inc * mult + inc) & _MASK128
+        inc = ((mult + 1) * inc) & _MASK128
+        mult = (mult * mult) & _MASK128
+        k >>= 1
+    return acc_mult, acc_inc
+
+
+def _pcg_advance(state, inc, k):
+    """The state of the LCG after `k` steps."""
+
+    mult, add = _pcg_power(inc, k & _MASK128)
+    return (mult * state + add) & _MASK128
+
+
+def _pcg_output(state):
+    """The output of PCG64DXSM for the state before a step."""
+
+    hi = state >> 64
+    lo = (state & _MASK64) | 1
+    hi ^= hi >> 32
+    hi = (hi * _PCG_MULT) & _MASK64
+    hi ^= hi >> 48
+    return (hi * lo) & _MASK64
+
+
+class PCG64DXSM(BitGenerator):
+    """
+    NumPy's PCG64DXSM bit generator on the device.
+
+    The permuted congruential generator of :class:`numpy.random.PCG64DXSM`,
+    a 128-bit LCG with the DXSM output function. It is seeded and has its
+    state as NumPy's, so the two can exchange states.
+
+    Parameters
+    ----------
+    seed : {None, int, array_like[ints], SeedSequence}, optional
+        The entropy the state is derived from, as by
+        :class:`numpy.random.PCG64DXSM`. See
+        :class:`dpnp.random.BitGenerator`.
+
+        Default: ``None``.
+    device : {None, string, SyclDevice, SyclQueue, Device}, optional
+        See :class:`dpnp.random.BitGenerator`.
+
+        Default: ``None``.
+    sycl_queue : {None, SyclQueue}, optional
+        See :class:`dpnp.random.BitGenerator`.
+
+        Default: ``None``.
+
+    Notes
+    -----
+    The raw outputs of :meth:`random_raw`, the floating-point values of
+    :meth:`dpnp.random.Generator.random` and
+    :meth:`dpnp.random.Generator.uniform`, and :meth:`dpnp.random.Generator.
+    bytes` take a fixed number of engine outputs for every value. They are
+    drawn in the order of the stream, and are the same as NumPy's for the
+    same state, where NumPy draws them with the same algorithm.
+
+    Every other method is sampled as by :meth:`jumped`: each work-item of a
+    call producing ``n`` values draws as many values as its rejection
+    samplers need from its own stream, the ``w``-th of which is that of
+    ``jumped(w)``, and the call advances the state as ``jumped(ceil(n /
+    8))``, but keeps the 32-bit half NumPy buffers. The values therefore
+    depend on the output size, and are not the same as NumPy's. For a given
+    state, output size and data type, the raw bits and the uniform floats
+    are the same on every device, while the samplers built on
+    transcendental functions agree to a few ULP, or differ where a
+    rejection sampler then takes another candidate.
+
+    """
+
+    _engine = ri.engine_id.pcg64dxsm
+
+    def __init__(self, seed=None, *, device=None, sycl_queue=None):
+        super().__init__(seed, device=device, sycl_queue=sycl_queue)
+
+        state = numpy.random.PCG64DXSM(self._seed_seq).state
+        self._state = state["state"]["state"]
+        self._inc = state["state"]["inc"]
+        self._has_uint32 = 0
+        self._uinteger = 0
+        self._jump = _pcg_power(self._inc, _PCG_JUMP)
+
+    def _set_state(self, state, inc, has_uint32, uinteger):
+        self._state = state
+        self._inc = inc
+        self._has_uint32 = has_uint32
+        self._uinteger = uinteger
+        self._jump = _pcg_power(inc, _PCG_JUMP)
+
+    @property
+    def state(self):
+        """
+        The state of the bit generator.
+
+        A dictionary ``{"bit_generator": "PCG64DXSM", "state": {"state":
+        state, "inc": inc}, "has_uint32": has_uint32, "uinteger":
+        uinteger}``, as that of :class:`numpy.random.PCG64DXSM`.
+
+        """
+
+        with self._lock:
+            return {
+                "bit_generator": self.__class__.__name__,
+                "state": {"state": self._state, "inc": self._inc},
+                "has_uint32": self._has_uint32,
+                "uinteger": self._uinteger,
+            }
+
+    @state.setter
+    def state(self, value):
+        if not isinstance(value, dict):
+            raise TypeError("state must be a dict")
+
+        bitgen = value.get("bit_generator", None)
+        if bitgen != self.__class__.__name__:
+            raise ValueError(
+                f"state must be for a {self.__class__.__name__} PRNG"
+            )
+
+        state = int(value["state"]["state"]) & _MASK128
+        inc = int(value["state"]["inc"]) & _MASK128
+        has_uint32 = int(value["has_uint32"])
+        uinteger = int(value["uinteger"]) & 0xFFFFFFFF
+        with self._lock:
+            self._set_state(state, inc, has_uint32, uinteger)
+
+    def advance(self, delta):
+        """
+        Advance the state as if `delta` raw values had been drawn.
+
+        For full documentation refer to
+        :obj:`numpy.random.PCG64DXSM.advance`.
+
+        """
+
+        delta = operator.index(delta)
+        with self._lock:
+            self._state = _pcg_advance(self._state, self._inc, delta)
+            self._has_uint32 = 0
+            self._uinteger = 0
+        return self
+
+    def jumped(self, jumps=1):
+        """
+        Return a new bit generator with the state jumped ahead.
+
+        For full documentation refer to
+        :obj:`numpy.random.PCG64DXSM.jumped`.
+
+        """
+
+        bit_generator = self.__class__(
+            self._seed_seq, sycl_queue=self._sycl_queue
+        )
+        bit_generator.state = self.state
+        bit_generator.advance(_PCG_JUMP * operator.index(jumps))
+        return bit_generator
+
+    def _words(self):
+        return [
+            self._state & _MASK64,
+            self._state >> 64,
+            self._inc & _MASK64,
+            self._inc >> 64,
+            self._jump[0] & _MASK64,
+            self._jump[0] >> 64,
+            self._jump[1] & _MASK64,
+            self._jump[1] >> 64,
+            self._has_uint32,
+            self._uinteger,
+        ]
+
+    def _launch(self, kernel, dtype, size):
+        halves = _PCG_SEQUENTIAL.get(kernel, {}).get(dtype, 0)
+        with self._lock:
+            words = self._words()
+            state, inc = self._state, self._inc
+            if halves == 2:
+                state = _pcg_advance(state, inc, size)
+            elif halves == 1:
+                if size and self._has_uint32:
+                    size -= 1
+                    self._has_uint32 = 0
+                raws = -(-size // 2)
+                if raws:
+                    # as NumPy, keep the high half of the last output, which
+                    # the next 32-bit draw takes if `size` is odd
+                    state = _pcg_advance(state, inc, raws - 1)
+                    self._uinteger = _pcg_output(state) >> 32
+                    self._has_uint32 = size % 2
+                    state = _pcg_advance(state, inc, 1)
+            else:
+                # work-item w draws from the stream jumped w times
+                n_wi = self._num_work_items(size)
+                state = _pcg_advance(state, inc, _PCG_JUMP * n_wi)
+            self._state = state
+        return words, 0, 0
+
+
+# L'Ecuyer's MRG32k3a: the moduli of its two components, and their
+# transition matrices, which map the state (x[n-3], x[n-2], x[n-1]) of a
+# component to the next one
+_MRG_M1 = 2**32 - 209
+_MRG_M2 = 2**32 - 22853
+_MRG_A1 = ((0, 1, 0), (0, 0, 1), (_MRG_M1 - 810728, 1403580, 0))
+_MRG_A2 = ((0, 1, 0), (0, 0, 1), (_MRG_M2 - 1370589, 0, 527612))
+
+# A substream is 2**76 values long
+_MRG_SUBSTREAM_BITS = 76
+
+
+def _mat_mul(a, b, m):
+    return tuple(
+        tuple(sum(a[i][k] * b[k][j] for k in range(3)) % m for j in range(3))
+        for i in range(3)
+    )
+
+
+def _mat_vec(a, v, m):
+    return tuple(sum(a[i][k] * v[k] for k in range(3)) % m for i in range(3))
+
+
+@functools.cache
+def _mrg_substream_jump(k):
+    """
+    The transition matrices of both components to the start of the
+    ``2**k``-th next substream.
+
+    """
+
+    if k == 0:
+        a1, a2 = _MRG_A1, _MRG_A2
+        for _ in range(_MRG_SUBSTREAM_BITS):
+            a1 = _mat_mul(a1, a1, _MRG_M1)
+            a2 = _mat_mul(a2, a2, _MRG_M2)
+        return a1, a2
+    a1, a2 = _mrg_substream_jump(k - 1)
+    return _mat_mul(a1, a1, _MRG_M1), _mat_mul(a2, a2, _MRG_M2)
+
+
+def _mrg_skip_substreams(s1, s2, n):
+    """The state `(s1, s2)` moved to the start of the `n`-th next substream."""
+
+    k = 0
+    while n:
+        if n & 1:
+            a1, a2 = _mrg_substream_jump(k)
+            s1 = _mat_vec(a1, s1, _MRG_M1)
+            s2 = _mat_vec(a2, s2, _MRG_M2)
+        n >>= 1
+        k += 1
+    return s1, s2
+
+
+class MRG32k3a(BitGenerator):
+    """
+    An MRG32k3a combined multiple recursive bit generator on the device.
+
+    The generator of L'Ecuyer (1999), with a period of about 2**191, as
+    implemented by the oneMKL device RNG API. Its state is two components
+    of three 32-bit words each.
+
+    Parameters
+    ----------
+    seed : {None, int, array_like[ints], SeedSequence}, optional
+        The entropy the state is derived from. See
+        :class:`dpnp.random.BitGenerator`.
+
+        Default: ``None``.
+    device : {None, string, SyclDevice, SyclQueue, Device}, optional
+        See :class:`dpnp.random.BitGenerator`.
+
+        Default: ``None``.
+    sycl_queue : {None, SyclQueue}, optional
+        See :class:`dpnp.random.BitGenerator`.
+
+        Default: ``None``.
+
+    Notes
+    -----
+    The engine gives values in ``[1, 2**32 - 209]``, which are spread over
+    the 32-bit range to give the raw 32-bit words, so that 209 words never
+    occur. A raw 64-bit value is two words, the first in the high half.
+
+    The stream is split into substreams of 2**76 values, as in L'Ecuyer et
+    al. (2002). Every call that produces ``n`` values takes ``ceil(n /
+    128)`` substreams, one for each work-item of the kernel, and every
+    work-item draws as many values as its rejection samplers need from its
+    own substream. The values therefore depend on the output size, and are
+    not the same as CuPy's. For a given state, output size and data type,
+    the raw bits and the uniform floats are the same on every device, while
+    the samplers built on transcendental functions agree to a few ULP, or
+    differ where a rejection sampler then takes another candidate.
+
+    References
+    ----------
+    .. [1] P. L'Ecuyer, "Good Parameters and Implementations for Combined
+           Multiple Recursive Random Number Generators", Operations
+           Research, 47(1), pp. 159-164, 1999.
+    .. [2] P. L'Ecuyer, R. Simard, E. J. Chen and W. D. Kelton, "An
+           Object-Oriented Random-Number Package with Many Long Streams and
+           Substreams", Operations Research, 50(6), pp. 1073-1075, 2002.
+
+    """
+
+    _engine = ri.engine_id.mrg32k3a
+
+    def __init__(self, seed=None, *, device=None, sycl_queue=None):
+        super().__init__(seed, device=device, sycl_queue=sycl_queue)
+
+        # values in [1, m - 1], which can't make a component all zero
+        words = [int(w) for w in self._seed_seq.generate_state(6, numpy.uint64)]
+        self._s1 = tuple(1 + w % (_MRG_M1 - 1) for w in words[:3])
+        self._s2 = tuple(1 + w % (_MRG_M2 - 1) for w in words[3:])
+
+    @staticmethod
+    def _check_component(value, m, name):
+        try:
+            words = tuple(operator.index(w) for w in value)
+        except TypeError:
+            raise TypeError(
+                f"{name} must be a sequence of 3 integers"
+            ) from None
+        if len(words) != 3:
+            raise ValueError(f"{name} must be a sequence of 3 integers")
+        if any(w < 0 or w >= m for w in words) or not any(words):
+            raise ValueError(
+                f"{name} must be in the range [0, {m}) and not all zero"
+            )
+        return words
+
+    @property
+    def state(self):
+        """
+        The state of the bit generator.
+
+        A dictionary ``{"bit_generator": "MRG32k3a", "state": {"s1": s1,
+        "s2": s2}}``, where `s1` and `s2` are the lists of the three words
+        of each component, the oldest first. The words of `s1` are less than
+        ``2**32 - 209``, those of `s2` less than ``2**32 - 22853``, and
+        neither is all zero.
+
+        """
+
+        with self._lock:
+            return {
+                "bit_generator": self.__class__.__name__,
+                "state": {"s1": list(self._s1), "s2": list(self._s2)},
+            }
+
+    @state.setter
+    def state(self, value):
+        if not isinstance(value, dict):
+            raise TypeError("state must be a dict")
+
+        bitgen = value.get("bit_generator", None)
+        if bitgen != self.__class__.__name__:
+            raise ValueError(
+                f"state must be for a {self.__class__.__name__} PRNG"
+            )
+
+        state = value["state"]
+        s1 = self._check_component(state["s1"], _MRG_M1, "s1")
+        s2 = self._check_component(state["s2"], _MRG_M2, "s2")
+        with self._lock:
+            self._s1 = s1
+            self._s2 = s2
+
+    def _launch(self, kernel, dtype, size):
+        # work-item w draws from the w-th next substream
+        n_wi = self._num_work_items(size)
+        with self._lock:
+            words = list(self._s1 + self._s2)
+            self._s1, self._s2 = _mrg_skip_substreams(self._s1, self._s2, n_wi)
+        return words, 0, 0
diff --git a/dpnp/random/_generator.py b/dpnp/random/_generator.py
index 89e690af643..a126870d5a5 100644
--- a/dpnp/random/_generator.py
+++ b/dpnp/random/_generator.py
@@ -332,7 +332,8 @@ class Generator:
     Parameters
     ----------
     bit_generator : BitGenerator
-        The bit generator, such as :class:`dpnp.random.Philox4x3210`.
+        The bit generator, such as :class:`dpnp.random.Philox4x3210` or
+        :class:`dpnp.random.PCG64DXSM`.
 
     Limitations
     -----------
@@ -344,8 +345,9 @@ class Generator:
     takes the USM type of the array parameters, or ``"device"`` when there
     are none.
 
-    Pickling, copying, and the ``advance`` and ``jumped`` methods of NumPy's
-    bit generators are not supported.
+    Pickling and copying are not supported, and only
+    :class:`dpnp.random.PCG64DXSM` has the ``advance`` and ``jumped``
+    methods of NumPy's bit generators.
 
     Parameters on the host (scalars, lists and NumPy arrays) are checked as
     in NumPy. Parameters on the device are not, since that would wait for
diff --git a/dpnp/tests/AGENTS.md b/dpnp/tests/AGENTS.md
index 81fa736ac91..640f7238b74 100644
--- a/dpnp/tests/AGENTS.md
+++ b/dpnp/tests/AGENTS.md
@@ -20,7 +20,7 @@
 - Error cases run on both libraries: parametrize over `xp` in `[dpnp, numpy]` so both must raise. Gate NumPy-version-dependent behaviour with `numpy_version()` (or `testing.with_requires` in CuPy tests).
 - Compute-follows-data error tests expect `(ValueError, ExecutionPlacementError)`: dpnp raises either, depending on the helper.
 - A new function is added to the matching parametrized lists in `test_sycl_queue.py` and `test_usm_type.py` (`test_1in_1out`, `test_2in_1out`, ...); elementwise ones also to `test_strides.py`.
-- `dpnp.random.Generator` values differ from NumPy's; test the samplers statistically, as the `test_moments` tests do. Cover the fp32 path of the discrete samplers with the `discrete_precision` fixture of `test_random_generator.py`, and restore any NumPy global random state a test changes, even when NumPy raises.
+- `dpnp.random` matches NumPy's values only for `PCG64DXSM` raw values, `random`, `uniform` and `bytes`; test the other samplers statistically, as the `test_moments` tests do. Cover the fp32 path of the discrete samplers with the `discrete_precision` fixture of `test_random_generator.py`, and restore any NumPy global random state a test changes, even when NumPy raises.
 - A test that may exceed the device's maximum allocation size skips on `USMAllocationError` (from `dpnp.exceptions`), not `MemoryError`, which it does not derive from.
 - Markers are strict and only `slow` and `multi_gpu` are registered. Tests that need huge allocations are `@pytest.mark.slow`; they never run in CI by design.
 - Skip lists are for device- or backend-specific failures only; use the full parametrized node ID relative to `dpnp/` (`tests/test_x.py::TestY::test_z[param]`; `conftest.py` matches it exactly, so a `dpnp/tests/...` entry matches nothing), and update or drop the entries of a test you rename or remove.
diff --git a/dpnp/tests/test_random_generator.py b/dpnp/tests/test_random_generator.py
index d8a680b93c2..262da7d5246 100644
--- a/dpnp/tests/test_random_generator.py
+++ b/dpnp/tests/test_random_generator.py
@@ -15,7 +15,13 @@
 
 import dpnp
 import dpnp.libtensor._rng_impl as ri
-from dpnp.random import Generator, Philox4x3210, default_rng
+from dpnp.random import (
+    PCG64DXSM,
+    Generator,
+    MRG32k3a,
+    Philox4x3210,
+    default_rng,
+)
 
 from .helper import get_float_dtypes, has_support_aspect64
 
@@ -33,8 +39,9 @@
 ]
 
 # the number of output elements each work-item of a kernel produces, with
-# Philox4x3210
+# Philox4x3210 and PCG64DXSM, and with MRG32k3a
 ITEMS_PER_WORK_ITEM = 8
+MRG_ITEMS_PER_WORK_ITEM = 128
 
 MASK32 = 0xFFFFFFFF
 
@@ -76,6 +83,66 @@ def _reference_raw(key, counter, n):
     return res
 
 
+MRG_M1 = 2**32 - 209
+MRG_M2 = 2**32 - 22853
+MRG_A1 = ((0, 1, 0), (0, 0, 1), (MRG_M1 - 810728, 1403580, 0))
+MRG_A2 = ((0, 1, 0), (0, 0, 1), (MRG_M2 - 1370589, 0, 527612))
+
+
+def _mrg32k3a_values(s1, s2):
+    """
+    Reference MRG32k3a (L'Ecuyer, 1999): the values in [1, m1] from the
+    state (s1, s2), as the oneMKL host API gives them.
+
+    """
+
+    x, y = list(s1), list(s2)
+    while True:
+        xn = (1403580 * x[1] - 810728 * x[0]) % MRG_M1
+        yn = (527612 * y[2] - 1370589 * y[0]) % MRG_M2
+        x = [x[1], x[2], xn]
+        y = [y[1], y[2], yn]
+        yield (xn - yn) % MRG_M1 or MRG_M1
+
+
+def _mrg_skip(a, v, m, k):
+    """The state `v` of a component with transition `a` after `k` steps."""
+
+    def mul(a, b):
+        return [
+            [sum(a[i][t] * b[t][j] for t in range(3)) % m for j in range(3)]
+            for i in range(3)
+        ]
+
+    p = [[int(i == j) for j in range(3)] for i in range(3)]
+    while k:
+        if k & 1:
+            p = mul(p, a)
+        a = mul(a, a)
+        k >>= 1
+    return [sum(p[i][t] * v[t] for t in range(3)) % m for i in range(3)]
+
+
+def _mrg_reference_raw(s1, s2, n):
+    """The expected output of MRG32k3a.random_raw(n) from state (s1, s2)."""
+
+    def word(r):
+        v = r - 1
+        return v + ((v * 209) >> 32)
+
+    n_wi = -(-n // MRG_ITEMS_PER_WORK_ITEM)
+    res = numpy.empty(n, dtype=numpy.uint64)
+    for w in range(n_wi):
+        k = w << 76
+        values = _mrg32k3a_values(
+            _mrg_skip(MRG_A1, s1, MRG_M1, k), _mrg_skip(MRG_A2, s2, MRG_M2, k)
+        )
+        for i in range(w, n, n_wi):
+            hi, lo = word(next(values)), word(next(values))
+            res[i] = (hi << 32) | lo
+    return res
+
+
 def _all_devices():
     return [
         d for d in dpctl.get_devices() if d.device_type.name in ("cpu", "gpu")
@@ -279,6 +346,279 @@ def test_param_overlaps_result(self):
         assert dpnp.isfinite(dpnp.asarray(x)).all()
 
 
+class TestPCG64DXSM:
+    @staticmethod
+    def _pair(seed):
+        bg = PCG64DXSM(seed)
+        ref = numpy.random.PCG64DXSM(seed)
+        assert bg.state == ref.state
+        return bg, ref
+
+    @pytest.mark.parametrize("seed", [0, 1234, (1 << 70) + 5])
+    @pytest.mark.parametrize("n", [0, 1, 7, 8, 9, 100, 1001])
+    def test_random_raw(self, seed, n):
+        bg, ref = self._pair(seed)
+        assert_array_equal(dpnp.asnumpy(bg.random_raw(n)), ref.random_raw(n))
+        assert bg.state == ref.state
+
+        bg.random_raw(n, output=False)
+        ref.random_raw(n, output=False)
+        assert bg.state == ref.state
+
+    @pytest.mark.parametrize("dtype", float_dtypes)
+    @pytest.mark.parametrize("n", [1, 5, 8, 17, 1000])
+    def test_random(self, dtype, n):
+        bg, ref = self._pair(3)
+        rng, ref_rng = Generator(bg), numpy.random.Generator(ref)
+        for _ in range(3):
+            res = rng.random(n, dtype=dtype)
+            assert_array_equal(dpnp.asnumpy(res), ref_rng.random(n, dtype))
+            assert bg.state == ref.state
+
+    @pytest.mark.skipif(not has_support_aspect64(), reason="fp64 required")
+    def test_uniform(self):
+        bg, ref = self._pair(4)
+        rng, ref_rng = Generator(bg), numpy.random.Generator(ref)
+        res = rng.uniform(-3.0, 7.5, 1001)
+        assert_array_equal(dpnp.asnumpy(res), ref_rng.uniform(-3.0, 7.5, 1001))
+
+        low = numpy.arange(10.0)
+        res = rng.uniform(dpnp.asarray(low), 20.0, (3, 10))
+        assert_array_equal(
+            dpnp.asnumpy(res), ref_rng.uniform(low, 20.0, (3, 10))
+        )
+        assert bg.state == ref.state
+
+    def test_bytes(self):
+        bg, ref = self._pair(5)
+        rng, ref_rng = Generator(bg), numpy.random.Generator(ref)
+        for length in [0, 1, 3, 4, 5, 13, -1, -3, -5, 1000]:
+            assert rng.bytes(length) == ref_rng.bytes(length)
+            assert bg.state == ref.state
+
+    @pytest.mark.parametrize("seed", range(5))
+    def test_mixed_calls(self, seed):
+        bg, ref = self._pair(seed)
+        rng, ref_rng = Generator(bg), numpy.random.Generator(ref)
+        sizes = numpy.random.default_rng(seed).integers(0, 40, size=20).tolist()
+        for i, n in enumerate(sizes):
+            kind = i % 4
+            if kind == 0:
+                res, exp = bg.random_raw(n), ref.random_raw(n)
+            elif kind == 1:
+                res = rng.random(n, dtype=dpnp.float32)
+                exp = ref_rng.random(n, dtype=numpy.float32)
+            elif kind == 2:
+                res, exp = rng.bytes(n), ref_rng.bytes(n)
+            else:
+                bg.advance(n)
+                ref.advance(n)
+                res = exp = None
+            if res is not None and not isinstance(res, bytes):
+                res = dpnp.asnumpy(res)
+            assert_array_equal(res, exp)
+            assert bg.state == ref.state
+
+    def test_rejection_sampler_streams(self):
+        # work-item w of a sampler whose draws vary draws from the stream of
+        # jumped(w), and the state advances as jumped(n_wi)
+        bg, ref = self._pair(6)
+        rng = Generator(bg)
+        rng.random(1, dtype=dpnp.float32)
+        state = bg.state
+
+        n_wi = 5
+        res = rng.integers(0, 1000, size=n_wi * ITEMS_PER_WORK_ITEM)
+        exp = numpy.empty(n_wi * ITEMS_PER_WORK_ITEM, dtype=numpy.int64)
+        for w in range(n_wi):
+            ref.state = state
+            ref_rng = numpy.random.Generator(ref.jumped(w))
+            exp[w::n_wi] = ref_rng.integers(0, 1000, ITEMS_PER_WORK_ITEM)
+        assert_array_equal(dpnp.asnumpy(res), exp)
+
+        ref.state = state
+        jumped = ref.jumped(n_wi).state
+        assert bg.state["state"] == jumped["state"]
+        assert bg.state["has_uint32"] == state["has_uint32"] == 1
+        assert bg.state["uinteger"] == state["uinteger"]
+
+    def test_state_exchange(self):
+        ref = numpy.random.PCG64DXSM(8)
+        numpy.random.Generator(ref).random(3, dtype=numpy.float32)
+        bg = PCG64DXSM()
+        bg.state = ref.state
+        assert bg.state == ref.state
+        assert_array_equal(dpnp.asnumpy(bg.random_raw(9)), ref.random_raw(9))
+
+        ref.state = bg.state
+        assert_array_equal(dpnp.asnumpy(bg.random_raw(9)), ref.random_raw(9))
+
+    def test_advance_jumped(self):
+        bg, ref = self._pair(9)
+        bg.random_raw(3)
+        ref.random_raw(3)
+        assert bg.advance(12345) is bg
+        ref.advance(12345)
+        assert bg.state == ref.state
+
+        jumped = bg.jumped(3)
+        assert jumped is not bg
+        assert jumped.state == ref.jumped(3).state
+        assert jumped.sycl_queue == bg.sycl_queue
+        assert bg.state == ref.state
+
+    def test_state_errors(self):
+        bg = PCG64DXSM(42)
+        assert_raises(TypeError, setattr, bg, "state", 1)
+        assert_raises(
+            ValueError, setattr, bg, "state", {"bit_generator": "PCG64"}
+        )
+
+    def test_spawn(self):
+        children = PCG64DXSM(5).spawn(2)
+        ref = numpy.random.PCG64DXSM(5).spawn(2)
+        for c, r in zip(children, ref):
+            assert c.state == r.state
+
+
+class TestMRG32k3a:
+    @staticmethod
+    def _with_state(s1, s2):
+        bg = MRG32k3a()
+        bg.state = {"bit_generator": "MRG32k3a", "state": {"s1": s1, "s2": s2}}
+        return bg
+
+    @pytest.mark.parametrize("n", [1, 127, 128, 129, 2000])
+    @pytest.mark.parametrize(
+        "s1, s2",
+        [
+            ([12345] * 3, [12345] * 3),
+            ([0, 0, 1], [MRG_M2 - 1, 0, 0]),
+            ([MRG_M1 - 1] * 3, [MRG_M2 - 1] * 3),
+        ],
+    )
+    def test_known_answer(self, n, s1, s2):
+        res = self._with_state(s1, s2).random_raw(n)
+        assert res.dtype == dpnp.uint64
+        assert_array_equal(dpnp.asnumpy(res), _mrg_reference_raw(s1, s2, n))
+
+    def test_state_advance(self):
+        bg = MRG32k3a(1234)
+        state = bg.state
+        assert state["bit_generator"] == "MRG32k3a"
+        s1, s2 = state["state"]["s1"], state["state"]["s2"]
+
+        # a draw of 257 values takes 3 substreams
+        bg.random_raw(257)
+        bg.random_raw(output=False)
+        k = 4 << 76
+        assert bg.state["state"] == {
+            "s1": _mrg_skip(MRG_A1, s1, MRG_M1, k),
+            "s2": _mrg_skip(MRG_A2, s2, MRG_M2, k),
+        }
+
+        # an empty draw does not advance the state
+        state = bg.state
+        bg.random_raw((0, 5))
+        assert bg.state == state
+
+    def test_state_roundtrip(self):
+        bg = MRG32k3a(42)
+        bg.random_raw(10)
+        state = bg.state
+        a = bg.random_raw(50)
+        bg.state = state
+        b = bg.random_raw(50)
+        assert_array_equal(a, b)
+
+    @pytest.mark.parametrize(
+        "s1, s2",
+        [
+            ([1, 2], [1, 2, 3]),
+            ([1, 2, 3], [1, 2, 3, 4]),
+            ([0, 0, 0], [1, 2, 3]),
+            ([1, 2, 3], [0, 0, 0]),
+            ([MRG_M1, 2, 3], [1, 2, 3]),
+            ([1, 2, 3], [MRG_M2, 2, 3]),
+            ([-1, 2, 3], [1, 2, 3]),
+        ],
+    )
+    def test_state_value_errors(self, s1, s2):
+        assert_raises(ValueError, self._with_state, s1, s2)
+
+    def test_state_errors(self):
+        bg = MRG32k3a(42)
+        assert_raises(TypeError, setattr, bg, "state", 1)
+        assert_raises(
+            ValueError, setattr, bg, "state", {"bit_generator": "Philox4x3210"}
+        )
+        assert_raises(TypeError, self._with_state, 1, [1, 2, 3])
+        assert_raises(TypeError, self._with_state, [1.5, 2, 3], [1, 2, 3])
+
+    def test_seed_sequence(self):
+        ss = numpy.random.SeedSequence(7)
+        bg = MRG32k3a(ss)
+        assert bg.seed_seq is ss
+        words = [int(w) for w in ss.generate_state(6, numpy.uint64)]
+        assert bg.state["state"] == {
+            "s1": [1 + w % (MRG_M1 - 1) for w in words[:3]],
+            "s2": [1 + w % (MRG_M2 - 1) for w in words[3:]],
+        }
+        assert MRG32k3a(7).state == bg.state
+
+    def test_spawn(self):
+        bg = MRG32k3a(5)
+        children = bg.spawn(2)
+        states = {str(c.state) for c in children}
+        states.add(str(bg.state))
+        assert len(states) == 3
+        for c in children:
+            assert c.sycl_queue == bg.sycl_queue
+
+    @pytest.mark.parametrize(
+        "name, params, mean, var",
+        [
+            ("random", (), 0.5, 1 / 12),
+            ("standard_normal", (), 0.0, 1.0),
+            ("standard_exponential", (), 1.0, 1.0),
+            ("standard_gamma", (0.5,), 0.5, 0.5),
+            ("integers", (0, 10), 4.5, 99 / 12),
+            ("poisson", (3.0,), 3.0, 3.0),
+            ("binomial", (20, 0.3), 6.0, 4.2),
+        ],
+    )
+    def test_moments(self, name, params, mean, var):
+        n = 10**6
+        res = getattr(Generator(MRG32k3a(3)), name)(*params, size=n)
+        res = dpnp.asnumpy(res).astype(numpy.float64)
+        assert abs(res.mean() - mean) <= 5 * numpy.sqrt(var / n)
+        assert abs(res.var() - var) <= 5 * var * numpy.sqrt(10 / n)
+
+
+@pytest.mark.parametrize("bitgen", [PCG64DXSM, MRG32k3a])
+class TestBitGeneratorState:
+    def test_lock_is_reentrant(self, bitgen):
+        bg = bitgen(1)
+
+        def _use_state():
+            with bg.lock:
+                bg.state = bg.state
+
+        t = threading.Thread(target=_use_state)
+        t.start()
+        t.join(timeout=10)
+        assert not t.is_alive()
+
+    @pytest.mark.parametrize("output", [True, False])
+    @pytest.mark.parametrize("size", [-8, (2, -8), 2.0])
+    def test_random_raw_bad_size(self, bitgen, size, output):
+        bg = bitgen(1)
+        bg.random_raw(9)
+        state = bg.state
+        assert_raises((TypeError, ValueError), bg.random_raw, size, output)
+        assert bg.state == state
+
+
 class TestDefaultRng:
     def test_reproducible(self):
         a = default_rng(123).random(100)
@@ -1869,7 +2209,7 @@ def _draw(self, device, bitgen):
             rng.hypergeometric(50, 70, 30, size=1000),
         )
 
-    @pytest.mark.parametrize("bitgen", [Philox4x3210])
+    @pytest.mark.parametrize("bitgen", [Philox4x3210, PCG64DXSM, MRG32k3a])
     def test_same_values(self, bitgen):
         devices = _all_devices()
         ref = [dpnp.asnumpy(x) for x in self._draw(devices[0], bitgen)]
diff --git a/dpnp/tests/test_sycl_queue.py b/dpnp/tests/test_sycl_queue.py
index 92bdeb52abc..a315c6ff6bb 100644
--- a/dpnp/tests/test_sycl_queue.py
+++ b/dpnp/tests/test_sycl_queue.py
@@ -988,7 +988,10 @@ def test_generator(func, args, kwargs, device):
     assert_sycl_queue_equal(res.sycl_queue, sycl_queue)
 
 
-@pytest.mark.parametrize("bit_generator", [dpnp.random.Philox4x3210])
+@pytest.mark.parametrize(
+    "bit_generator",
+    [dpnp.random.Philox4x3210, dpnp.random.PCG64DXSM, dpnp.random.MRG32k3a],
+)
 @pytest.mark.parametrize("device", valid_dev, ids=dev_ids)
 def test_bit_generator(bit_generator, device):
     sycl_queue = dpctl.SyclQueue(device, property="in_order")
diff --git a/dpnp/tests/test_usm_type.py b/dpnp/tests/test_usm_type.py
index 08577bde648..f9fccdc4072 100644
--- a/dpnp/tests/test_usm_type.py
+++ b/dpnp/tests/test_usm_type.py
@@ -1914,7 +1914,10 @@ def test_generator_permutations(usm_type):
     assert x.usm_type == usm_type
 
 
-@pytest.mark.parametrize("bit_generator", [dpnp.random.Philox4x3210])
+@pytest.mark.parametrize(
+    "bit_generator",
+    [dpnp.random.Philox4x3210, dpnp.random.PCG64DXSM, dpnp.random.MRG32k3a],
+)
 @pytest.mark.parametrize("usm_type", list_of_usm_types)
 def test_bit_generator(bit_generator, usm_type):
     res = bit_generator(1).random_raw(5, usm_type=usm_type)
diff --git a/dpnp/tests/third_party/cupy/random_tests/test_bit_generator.py b/dpnp/tests/third_party/cupy/random_tests/test_bit_generator.py
index 03b64074188..d5f5d1bdeaa 100644
--- a/dpnp/tests/third_party/cupy/random_tests/test_bit_generator.py
+++ b/dpnp/tests/third_party/cupy/random_tests/test_bit_generator.py
@@ -56,7 +56,6 @@ def setUp(self):
 
 @testing.with_requires("numpy>=1.17.0")
 @testing.fix_random()
-@pytest.mark.skip("MRG32k3a is not supported yet")
 class TestBitGeneratorMRG32k3a(BitGeneratorTestCase, unittest.TestCase):
     def setUp(self):
         super().setUp()