From e06fd8c3ca9a8d9b2197895dfc256209e022f8c7 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                      |  14 +-
 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                 | 403 ++++++++++++++++++
 dpnp/random/_generator.py                     |   3 +-
 dpnp/tests/AGENTS.md                          |   2 +-
 dpnp/tests/test_random_generator.py           | 322 +++++++++++++-
 dpnp/tests/test_sycl_queue.py                 |   5 +-
 dpnp/tests/test_usm_type.py                   |   5 +-
 16 files changed, 1082 insertions(+), 27 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1a8c8e024bc0..d7df4c2ab55f 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 b8f0ec87e530..64e2bb159d18 100644
--- a/doc/known_words.txt
+++ b/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
diff --git a/doc/reference/random.rst b/doc/reference/random.rst
index 0a7f4a300f39..810a5ff7f5bd 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.
@@ -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
diff --git a/dpnp/AGENTS.md b/dpnp/AGENTS.md
index f87cde807604..8d2a8d977af2 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 abdd80042eb1..4c4396807480 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 c9af142bd7d5..bc7791e9226e 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 481e461f79ce..ddaefdea3031 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 values 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.
@@ -53,8 +63,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 &&
@@ -63,9 +78,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;
 }
 
@@ -76,6 +93,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>);
 
@@ -103,7 +122,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 7e772e5e1dce..b98c0e90dc22 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 512243045ecb..0e85998bc486 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 4e7a50ad44d9..dbf860948850 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 ee90b2d6fb8d..b705fc6e9f01 100644
--- a/dpnp/random/_bit_generator.py
+++ b/dpnp/random/_bit_generator.py
@@ -40,6 +40,8 @@
 # pylint: disable=no-name-in-module
 # pylint: disable=protected-access
 
+import functools
+import operator
 import threading
 
 import dpctl.utils as dpu
@@ -50,6 +52,8 @@
 
 __all__ = [
     "BitGenerator",
+    "MRG32k3a",
+    "PCG64DXSM",
     "Philox4x3210",
     "SeedSequence",
 ]
@@ -391,3 +395,402 @@ def _launch(self, kernel, dtype, size):
             counter = self._counter
             self._counter = (counter + (n_blocks << 64)) & _MASK128
         return [self._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 only on the state and the output size, and are the same on every
+    device, but not the same as NumPy's.
+
+    """
+
+    _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 only on the state and the
+    output size, and are the same on every device, but not the same as
+    CuPy's.
+
+    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 af9fc8ef9344..4978d5ad53a9 100644
--- a/dpnp/random/_generator.py
+++ b/dpnp/random/_generator.py
@@ -301,7 +301,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
     -----------
diff --git a/dpnp/tests/AGENTS.md b/dpnp/tests/AGENTS.md
index 81fa736ac919..640f7238b74b 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 3a96e07f2f62..7212e2d120ab 100644
--- a/dpnp/tests/test_random_generator.py
+++ b/dpnp/tests/test_random_generator.py
@@ -12,7 +12,13 @@
 )
 
 import dpnp
-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
 
@@ -30,8 +36,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
 
@@ -73,6 +80,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")
@@ -158,6 +225,255 @@ def test_base_class(self):
         assert_raises(TypeError, dpnp.random.BitGenerator)
 
 
+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)
+
+
 class TestDefaultRng:
     def test_reproducible(self):
         a = default_rng(123).random(100)
@@ -1591,7 +1907,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 92bdeb52abc5..a315c6ff6bbe 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 fd48df0945cb..6550d1f2e3dc 100644
--- a/dpnp/tests/test_usm_type.py
+++ b/dpnp/tests/test_usm_type.py
@@ -1899,7 +1899,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)