diff --git a/doc/Jamfile.v2 b/doc/Jamfile.v2
index 7a7bc91c5..ba9e645bb 100644
--- a/doc/Jamfile.v2
+++ b/doc/Jamfile.v2
@@ -71,6 +71,7 @@ doxygen_files =
     traits
     splitmix64
     xoshiro
+    mwcx
     inverse_gaussian_distribution
     generalized_inverse_gaussian_distribution
 ;
diff --git a/doc/generators.qbk b/doc/generators.qbk
index f43b6ff6d..863a6741e 100644
--- a/doc/generators.qbk
+++ b/doc/generators.qbk
@@ -93,6 +93,10 @@ numbers mean faster random number generation.
   [[__xoshiro128pp] [2[sup 128]] [`4*sizeof(uint32_t)`] [133%] [xoshiro128++ from https://prng.di.unimi.it. Returns uint32_t]]
   [[__xoshiro128f] [2[sup 128]] [`4*sizeof(uint32_t)`] [145%] [This generator returns float instead of uint32_t. It is modified xoshiro128+ from https://prng.di.unimi.it]]
   [[__xoshiro128mm] [2[sup 128]] [`4*sizeof(uint32_t)`] [116%] [xoshiro128** from https://prng.di.unimi.it. Returns uint32_t]]
+  [[__mwc192x] [~2[sup 191]] [`3*sizeof(uint64_t)`] [-] [MWC192X, a scrambled multiply-with-carry generator from https://prng.di.unimi.it (see https://arxiv.org/abs/2610.27313). Besides discard(), provides jump(c, e), which advances the generator by c*2[sup e] steps, to obtain non-overlapping subsequences at distances larger than discard() can express: jump(1, 96) and jump(1, 144) are equivalent to the jump() and long_jump() functions of the reference code.]]
+  [[__mwc128x] [~2[sup 127]] [`2*sizeof(uint64_t)`] [-] [MWC128X, a scrambled multiply-with-carry generator from https://prng.di.unimi.it (see https://arxiv.org/abs/2610.27313). Besides discard(), provides jump(c, e), which advances the generator by c*2[sup e] steps, to obtain non-overlapping subsequences at distances larger than discard() can express: jump(1, 64) and jump(1, 96) are equivalent to the jump() and long_jump() functions of the reference code.]]
+  [[__mwc128x32] [~2[sup 127]] [`4*sizeof(uint32_t)`] [-] [MWC128X/32, a scrambled multiply-with-carry generator from https://prng.di.unimi.it (see https://arxiv.org/abs/2610.27313). Returns uint32_t. Besides discard(), provides jump(c, e), which advances the generator by c*2[sup e] steps, to obtain non-overlapping subsequences at distances larger than discard() can express: jump(1, 64) and jump(1, 96) are equivalent to the jump() and long_jump() functions of the reference code.]]
+  [[__mwc64x32] [~2[sup 63]] [`2*sizeof(uint32_t)`] [-] [MWC64X/32, a scrambled multiply-with-carry generator from https://prng.di.unimi.it (see https://arxiv.org/abs/2610.27313). Returns uint32_t. Besides discard(), provides jump(c, e), which advances the generator by c*2[sup e] steps, to obtain non-overlapping subsequences at distances larger than discard() can express: jump(1, 32) and jump(1, 48) are equivalent to the jump() and long_jump() functions of the reference code.]]
 
 ]
 
diff --git a/doc/random.qbk b/doc/random.qbk
index cbc5f91ec..09640aa54 100644
--- a/doc/random.qbk
+++ b/doc/random.qbk
@@ -79,6 +79,10 @@
 [def __xoshiro128pp [classref boost::random::xoshiro128pp xoshiro128pp]]
 [def __xoshiro128f [classref boost::random::xoshiro128f xoshiro128f]]
 [def __xoshiro128mm [classref boost::random::xoshiro128mm xoshiro128mm]]
+[def __mwc192x [classref boost::random::mwc192x mwc192x]]
+[def __mwc128x [classref boost::random::mwc128x mwc128x]]
+[def __mwc128x32 [classref boost::random::mwc128x32 mwc128x32]]
+[def __mwc64x32 [classref boost::random::mwc64x32 mwc64x32]]
 
 [def __uniform_smallint [classref boost::random::uniform_smallint uniform_smallint]]
 [def __uniform_int_distribution [classref boost::random::uniform_int_distribution uniform_int_distribution]]
diff --git a/include/boost/random.hpp b/include/boost/random.hpp
index 884aa93ec..695f4202d 100644
--- a/include/boost/random.hpp
+++ b/include/boost/random.hpp
@@ -50,6 +50,7 @@
 #include <boost/random/xor_combine.hpp>
 #include <boost/random/splitmix64.hpp>
 #include <boost/random/xoshiro.hpp>
+#include <boost/random/mwcx.hpp>
 
 // misc
 #include <boost/random/generate_canonical.hpp>
diff --git a/include/boost/random/detail/mwcx_arithmetic.hpp b/include/boost/random/detail/mwcx_arithmetic.hpp
new file mode 100644
index 000000000..5aac911b7
--- /dev/null
+++ b/include/boost/random/detail/mwcx_arithmetic.hpp
@@ -0,0 +1,211 @@
+/*
+ * Copyright Sebastiano Vigna 2026.
+ * Distributed under the Boost Software License, Version 1.0. (See
+ * accompanying file LICENSE_1_0.txt or copy at
+ * http://www.boost.org/LICENSE_1_0.txt)
+ *
+ * See http://www.boost.org for most recent version including documentation.
+ */
+
+#ifndef BOOST_RANDOM_DETAIL_MWCX_ARITHMETIC_HPP
+#define BOOST_RANDOM_DETAIL_MWCX_ARITHMETIC_HPP
+
+#include <boost/config.hpp>
+#include <array>
+#include <cstddef>
+#include <cstdint>
+
+// The step of the generators with 64-bit digits needs the 128-bit result of
+// a * x + c, where a, x, and c are 64-bit words. Defining
+// BOOST_RANDOM_MWCX_PORTABLE_MUL forces the portable implementation, which
+// is used when the compiler provides neither a 128-bit integer type nor a
+// 64 x 64 -> 128-bit multiplication intrinsic (e.g., on 32-bit targets).
+#if !defined(BOOST_RANDOM_MWCX_PORTABLE_MUL)
+#  if defined(__SIZEOF_INT128__)
+#    define BOOST_RANDOM_MWCX_HAS_INT128
+#  elif defined(_MSC_VER) && (defined(_M_X64) || defined(_M_AMD64)) && !defined(_M_ARM64EC)
+#    include <intrin.h>
+#    pragma intrinsic(_umul128)
+#    define BOOST_RANDOM_MWCX_HAS_UMUL128
+#  elif defined(_MSC_VER) && (defined(_M_ARM64) || defined(_M_ARM64EC))
+#    include <intrin.h>
+#    define BOOST_RANDOM_MWCX_HAS_UMULH
+#  endif
+#endif
+
+namespace boost {
+namespace random {
+namespace detail {
+
+/**
+ * Returns the lower 64 bits of a * x + c, and stores the upper 64 bits in hi.
+ *
+ * The sum cannot overflow 128 bits, as (2^64 - 1)^2 + 2^64 - 1 < 2^128.
+ */
+inline std::uint64_t mwcx_mul_add(const std::uint64_t a, const std::uint64_t x, const std::uint64_t c, std::uint64_t& hi) noexcept
+{
+#if defined(BOOST_RANDOM_MWCX_HAS_INT128)
+    __extension__ typedef unsigned __int128 uint128_type;
+    const uint128_type t = static_cast<uint128_type>(a) * x + c;
+    hi = static_cast<std::uint64_t>(t >> 64);
+    return static_cast<std::uint64_t>(t);
+#elif defined(BOOST_RANDOM_MWCX_HAS_UMUL128)
+    std::uint64_t h;
+    const std::uint64_t l = _umul128(a, x, &h) + c;
+    hi = h + (l < c);
+    return l;
+#elif defined(BOOST_RANDOM_MWCX_HAS_UMULH)
+    const std::uint64_t l = a * x + c;
+    hi = __umulh(a, x) + (l < c);
+    return l;
+#else
+    // Schoolbook multiplication on 32-bit halves. The two halves of c are
+    // added to the first two partial products: each partial product can
+    // absorb two 32-bit addends without overflow, as
+    // (2^32 - 1)^2 + 2 * (2^32 - 1) = 2^64 - 1.
+    const std::uint64_t mask = UINT64_C(0xffffffff);
+    const std::uint64_t a0 = a & mask, a1 = a >> 32, x0 = x & mask, x1 = x >> 32;
+    const std::uint64_t t0 = a0 * x0 + (c & mask);
+    const std::uint64_t t1 = a0 * x1 + (t0 >> 32) + (c >> 32);
+    const std::uint64_t t2 = a1 * x0 + (t1 & mask);
+    hi = a1 * x1 + (t1 >> 32) + (t2 >> 32);
+    return (t2 << 32) | (t0 & mask);
+#endif
+}
+
+// Minimal modular arithmetic on multiple-precision natural numbers, used by
+// discard() and jump(). Numbers are little-endian arrays of 64-bit words, and
+// operands are assumed to be smaller than the modulus. Multiplication is
+// shift-and-add, which needs no wide multiplication: it is slow, but jumps
+// are rare, and at these sizes a proper multiple-precision library would not
+// be much faster.
+
+/** Returns whether a >= b. */
+template <std::size_t N>
+inline bool mwcx_ge(const std::array<std::uint64_t, N>& a, const std::array<std::uint64_t, N>& b) noexcept
+{
+    for (std::size_t i = N; i-- != 0;)
+    {
+        if (a[i] != b[i])
+        {
+            return a[i] > b[i];
+        }
+    }
+
+    return true;
+}
+
+/** Returns whether 0 < s < m, that is, whether s is a valid state of an MWC generator with modulus m. */
+template <std::size_t N>
+inline bool mwcx_is_valid(const std::array<std::uint64_t, N>& s, const std::array<std::uint64_t, N>& m) noexcept
+{
+    bool nonzero = false;
+    for (std::size_t i = 0; i < N; ++i)
+    {
+        nonzero |= s[i] != 0;
+    }
+
+    return nonzero && !mwcx_ge(s, m);
+}
+
+/** Returns the number of bits of the little-endian integer of n words held in x (zero for zero). */
+inline std::size_t mwcx_bit_len(const std::uint64_t* x, std::size_t n) noexcept
+{
+    while (n != 0 && x[n - 1] == 0)
+    {
+        --n;
+    }
+
+    if (n == 0)
+    {
+        return 0;
+    }
+
+    std::size_t bits = 64 * n;
+    for (std::uint64_t w = x[n - 1]; (w & (UINT64_C(1) << 63)) == 0; w <<= 1)
+    {
+        --bits;
+    }
+
+    return bits;
+}
+
+/** Sets a to (a + b) mod m. */
+template <std::size_t N>
+inline void mwcx_add_mod(std::array<std::uint64_t, N>& a, const std::array<std::uint64_t, N>& b, const std::array<std::uint64_t, N>& m) noexcept
+{
+    std::uint64_t carry = 0;
+    for (std::size_t i = 0; i < N; ++i)
+    {
+        const std::uint64_t s = a[i] + b[i];
+        const std::uint64_t c0 = s < a[i];
+        a[i] = s + carry;
+        carry = c0 | (a[i] < s);
+    }
+
+    // The sum is smaller than 2m, so one subtraction suffices; if the sum
+    // overflowed, the subtraction wraps around to the correct result.
+    if (carry != 0 || mwcx_ge(a, m))
+    {
+        std::uint64_t borrow = 0;
+        for (std::size_t i = 0; i < N; ++i)
+        {
+            const std::uint64_t d = a[i] - m[i];
+            const std::uint64_t b0 = d > a[i];
+            a[i] = d - borrow;
+            borrow = b0 | (a[i] > d);
+        }
+    }
+}
+
+/** Returns a * b mod m. */
+template <std::size_t N>
+inline std::array<std::uint64_t, N> mwcx_mul_mod(const std::array<std::uint64_t, N>& a, const std::array<std::uint64_t, N>& b, const std::array<std::uint64_t, N>& m) noexcept
+{
+    std::array<std::uint64_t, N> r {};
+    std::array<std::uint64_t, N> t = a;
+    const std::size_t bits = mwcx_bit_len(b.data(), N);
+    for (std::size_t i = 0; i < bits; ++i)
+    {
+        if ((b[i / 64] >> (i % 64)) & 1)
+        {
+            mwcx_add_mod(r, t, m);
+        }
+
+        const std::array<std::uint64_t, N> u = t;
+        mwcx_add_mod(t, u, m);
+    }
+
+    return r;
+}
+
+/** Returns base^n mod m, where n is the little-endian integer of len words held in exp. */
+template <std::size_t N>
+inline std::array<std::uint64_t, N> mwcx_pow_mod(const std::array<std::uint64_t, N>& base, const std::uint64_t* exp, const std::size_t len, const std::array<std::uint64_t, N>& m) noexcept
+{
+    std::array<std::uint64_t, N> r {};
+    r[0] = 1;
+    std::array<std::uint64_t, N> b = base;
+    // Right-to-left square-and-multiply
+    const std::size_t bits = mwcx_bit_len(exp, len);
+    for (std::size_t i = 0; i < bits; ++i)
+    {
+        if ((exp[i / 64] >> (i % 64)) & 1)
+        {
+            r = mwcx_mul_mod(r, b, m);
+        }
+
+        if (i + 1 < bits)
+        {
+            b = mwcx_mul_mod(b, b, m);
+        }
+    }
+
+    return r;
+}
+
+} // namespace detail
+} // namespace random
+} // namespace boost
+
+#endif // BOOST_RANDOM_DETAIL_MWCX_ARITHMETIC_HPP
diff --git a/include/boost/random/detail/mwcx_base.hpp b/include/boost/random/detail/mwcx_base.hpp
new file mode 100644
index 000000000..b5207d575
--- /dev/null
+++ b/include/boost/random/detail/mwcx_base.hpp
@@ -0,0 +1,395 @@
+/*
+ * Copyright Sebastiano Vigna 2026.
+ * Distributed under the Boost Software License, Version 1.0. (See
+ * accompanying file LICENSE_1_0.txt or copy at
+ * http://www.boost.org/LICENSE_1_0.txt)
+ *
+ * See http://www.boost.org for most recent version including documentation.
+ */
+
+#ifndef BOOST_RANDOM_DETAIL_MWCX_BASE_HPP
+#define BOOST_RANDOM_DETAIL_MWCX_BASE_HPP
+
+#include <boost/random/detail/mwcx_arithmetic.hpp>
+#include <boost/throw_exception.hpp>
+#include <boost/config.hpp>
+#include <array>
+#include <cstddef>
+#include <cstdint>
+#include <ios>
+#include <istream>
+#include <iterator>
+#include <limits>
+#include <stdexcept>
+#include <type_traits>
+
+namespace boost {
+namespace random {
+namespace detail {
+
+/**
+ * Common base of the scrambled multiply-with-carry generators.
+ *
+ * Digit is the type of the digits (std::uint32_t or std::uint64_t), K is the
+ * number of words of the state (the digits and the carry), and N is the
+ * number of 64-bit words of the state. The state is stored in state_, least
+ * significant word first: the oldest digit is state_[0], and the carry is
+ * state_[K - 1]. Read as an integer s, the state is an element of the
+ * multiplicative group of the integers modulo the prime m = a * b^(K - 1) - 1,
+ * where b = 2^w is the base and w is the size in bits of a digit; each step
+ * multiplies s by b^(-1) mod m.
+ *
+ * Derived must provide next() and the static functions modulus() (m) and
+ * inverse_base() (b^(-1) mod m), returning std::array<std::uint64_t, N>.
+ */
+template <typename Derived, typename Digit, std::size_t K, std::size_t N>
+class mwcx_base
+{
+    static_assert(std::is_same<Digit, std::uint32_t>::value || std::is_same<Digit, std::uint64_t>::value,
+                  "The digits must be 32-bit or 64-bit unsigned integers");
+    static_assert(K * sizeof(Digit) == N * sizeof(std::uint64_t), "Inconsistent state size");
+
+protected:
+
+    using words_type = std::array<std::uint64_t, N>;
+
+    std::array<Digit, K> state_;
+
+private:
+
+    // Below this distance, discard() steps the generator instead of jumping.
+    // The cost of a jump grows quadratically with the size of the state, and
+    // linearly with the logarithm of the distance; this threshold is close to
+    // the break-even point on current x86-64 processors.
+    static constexpr std::uint64_t discard_threshold = UINT64_C(1024) << (2 * N);
+
+    static words_type to_words(const std::array<Digit, K>& s, const std::true_type&) noexcept
+    {
+        words_type w;
+        for (std::size_t i = 0; i < N; ++i)
+        {
+            w[i] = s[i];
+        }
+
+        return w;
+    }
+
+    static words_type to_words(const std::array<Digit, K>& s, const std::false_type&) noexcept
+    {
+        words_type w;
+        for (std::size_t i = 0; i < N; ++i)
+        {
+            w[i] = static_cast<std::uint64_t>(s[2 * i]) | static_cast<std::uint64_t>(s[2 * i + 1]) << 32;
+        }
+
+        return w;
+    }
+
+    static std::array<Digit, K> from_words(const words_type& w, const std::true_type&) noexcept
+    {
+        std::array<Digit, K> s;
+        for (std::size_t i = 0; i < K; ++i)
+        {
+            s[i] = static_cast<Digit>(w[i]);
+        }
+
+        return s;
+    }
+
+    static std::array<Digit, K> from_words(const words_type& w, const std::false_type&) noexcept
+    {
+        std::array<Digit, K> s;
+        for (std::size_t i = 0; i < N; ++i)
+        {
+            s[2 * i] = static_cast<Digit>(w[i]);
+            s[2 * i + 1] = static_cast<Digit>(w[i] >> 32);
+        }
+
+        return s;
+    }
+
+    using is_64 = std::integral_constant<bool, sizeof(Digit) == sizeof(std::uint64_t)>;
+
+    words_type words() const noexcept
+    {
+        return to_words(state_, is_64());
+    }
+
+    void set_words(const words_type& w) noexcept
+    {
+        state_ = from_words(w, is_64());
+    }
+
+    // Sets the state to s; if s is not a valid state, the carry is set to one.
+    void set_state(const words_type& w) noexcept
+    {
+        set_words(w);
+        if (!mwcx_is_valid(w, Derived::modulus()))
+        {
+            state_[K - 1] = 1;
+        }
+    }
+
+    // Multiplies the state, read as an integer, by r modulo m.
+    void jump_mul(const words_type& r) noexcept
+    {
+        set_words(mwcx_mul_mod(words(), r, Derived::modulus()));
+    }
+
+    static std::uint64_t splitmix(std::uint64_t& x) noexcept
+    {
+        std::uint64_t z = (x += UINT64_C(0x9E3779B97F4A7C15));
+        z = (z ^ (z >> 30)) * UINT64_C(0xBF58476D1CE4E5B9);
+        z = (z ^ (z >> 27)) * UINT64_C(0x94D049BB133111EB);
+        return z ^ (z >> 31);
+    }
+
+public:
+
+    using result_type = Digit;
+    using seed_type = std::uint64_t;
+
+    // Required for old Boost.Random concept
+    static constexpr bool has_fixed_range {false};
+
+    /**
+     * Seeds the generator with a user provided seed, as suggested by the
+     * reference implementation: the carry is set to one, and the digits are
+     * filled with the output of a SplitMix generator seeded with @c value
+     * (32-bit digits take the lower and then the upper half of each output).
+     */
+    void seed(const seed_type value) noexcept
+    {
+        std::uint64_t x = value;
+        std::uint64_t z = 0;
+        for (std::size_t i = 0; i < K - 1; ++i)
+        {
+            BOOST_IF_CONSTEXPR (sizeof(Digit) == sizeof(std::uint64_t))
+            {
+                state_[i] = static_cast<Digit>(splitmix(x));
+            }
+            else
+            {
+                if (i % 2 == 0)
+                {
+                    z = splitmix(x);
+                }
+
+                state_[i] = static_cast<Digit>(z >> (32 * (i % 2)));
+            }
+        }
+
+        state_[K - 1] = 1;
+    }
+
+    /** Seeds the generator with the default seed (zero). */
+    void seed() noexcept
+    {
+        seed(static_cast<seed_type>(0));
+    }
+
+    /**
+     * Sets the state of the generator to the little-endian integer formed by
+     * 32-bit values produced by @c seq.generate(). If the integer is not a
+     * valid state, that is, if it is zero or not smaller than the modulus,
+     * the carry is set to one.
+     */
+    template <typename Sseq, typename std::enable_if<!std::is_convertible<Sseq, seed_type>::value, bool>::type = true>
+    void seed(Sseq& seq)
+    {
+        std::array<std::uint32_t, 2 * N> seeds;
+        seq.generate(seeds.begin(), seeds.end());
+
+        words_type w;
+        for (std::size_t i = 0; i < N; ++i)
+        {
+            w[i] = static_cast<std::uint64_t>(seeds[2 * i]) | static_cast<std::uint64_t>(seeds[2 * i + 1]) << 32;
+        }
+
+        set_state(w);
+    }
+
+    /**
+     * Sets the state of the generator using values from an iterator range:
+     * the digits, oldest first, and then the carry. If the resulting state is
+     * not valid, that is, if it is zero or not smaller than the modulus, the
+     * carry is set to one.
+     */
+    template <typename FIter>
+    void seed(FIter first, FIter last)
+    {
+        static_assert(std::is_integral<typename std::iterator_traits<FIter>::value_type>::value,
+                      "Value type must be a built-in integer type" );
+
+        std::array<Digit, K> s;
+        std::size_t offset = 0;
+        while (first != last && offset < K)
+        {
+            s[offset++] = static_cast<Digit>(*first++);
+        }
+
+        if (offset != K)
+        {
+            boost::throw_exception(std::invalid_argument("Not enough elements in call to seed."));
+        }
+
+        set_state(to_words(s, is_64()));
+    }
+
+    /** Constructs a generator and calls @c seed(). */
+    mwcx_base() noexcept { seed(); }
+
+    /** Seeds the generator with a user provided seed. */
+    explicit mwcx_base(const seed_type value) noexcept
+    {
+        seed(value);
+    }
+
+    /** Sets the state of the generator using values from an iterator range. */
+    template <typename FIter>
+    mwcx_base(FIter& first, FIter last) { seed(first, last); }
+
+    /**
+     * Sets the state of the generator using 32-bit values produced by @c seq.generate().
+     *
+     * @xmlnote
+     * The copy constructor will always be preferred over
+     * the templated constructor.
+     * @endxmlnote
+     */
+    template <typename Sseq, typename std::enable_if<!std::is_convertible<Sseq, mwcx_base>::value, bool>::type = true>
+    explicit mwcx_base(Sseq& seq)
+    {
+        seed(seq);
+    }
+
+    inline result_type next() noexcept
+    {
+        return static_cast<Derived*>(this)->next();
+    }
+
+    /** Returns the next value of the generator. */
+    inline result_type operator()() noexcept
+    {
+        return next();
+    }
+
+    /**
+     * Advances the state of the generator by @c z steps. Large distances are
+     * computed in time logarithmic in @c z, by multiplying the state by
+     * b^(-z) mod m.
+     */
+    inline void discard(const std::uint64_t z) noexcept
+    {
+        if (z < discard_threshold)
+        {
+            for (std::uint64_t i = 0; i < z; ++i)
+            {
+                next();
+            }
+        }
+        else
+        {
+            jump(z, 0);
+        }
+    }
+
+    /**
+     * Advances the state of the generator by c * 2^e steps, in time
+     * logarithmic in c plus linear in e, by multiplying the state by
+     * b^(-c * 2^e) mod m. This makes it possible to obtain non-overlapping
+     * subsequences at distances larger than discard() can express. For the
+     * jump to be meaningful, c * 2^e should be smaller than the period of the
+     * generator.
+     */
+    inline void jump(const std::uint64_t c, const unsigned int e) noexcept
+    {
+        words_type r = mwcx_pow_mod(Derived::inverse_base(), &c, 1, Derived::modulus());
+        for (unsigned int i = 0; i < e; ++i)
+        {
+            r = mwcx_mul_mod(r, r, Derived::modulus());
+        }
+
+        jump_mul(r);
+    }
+
+    /**
+     * Returns true if the two generators will produce identical
+     * sequences of values.
+     */
+    inline friend bool operator==(const mwcx_base& lhs, const mwcx_base& rhs) noexcept
+    {
+        return lhs.state_ == rhs.state_;
+    }
+
+    /**
+     * Returns true if the two generators will produce different
+     * sequences of values.
+     */
+    inline friend bool operator!=(const mwcx_base& lhs, const mwcx_base& rhs) noexcept
+    {
+        return lhs.state_ != rhs.state_;
+    }
+
+    /** Writes the state of the generator to a @c std::ostream. */
+    template <typename CharT, typename Traits>
+    inline friend std::basic_ostream<CharT, Traits>& operator<<(std::basic_ostream<CharT, Traits>& ost,
+                                                                const mwcx_base& e)
+    {
+        for (std::size_t i = 0; i < K; ++i)
+        {
+            ost << e.state_[i] << ' ';
+        }
+
+        return ost;
+    }
+
+    /** Reads the state of the generator from a @c std::istream. */
+    template <typename CharT, typename Traits>
+    inline friend std::basic_istream<CharT, Traits>& operator>>(std::basic_istream<CharT, Traits>& ist,
+                                                                mwcx_base& e)
+    {
+        for (std::size_t i = 0; i < K; ++i)
+        {
+            ist >> e.state_[i] >> std::ws;
+        }
+
+        return ist;
+    }
+
+    /** Fills a range with random values */
+    template <typename FIter>
+    inline void generate(FIter first, FIter last) noexcept
+    {
+        using iter_type = typename std::iterator_traits<FIter>::value_type;
+
+        while (first != last)
+        {
+            *first++ = static_cast<iter_type>(next());
+        }
+    }
+
+    /** Returns the largest value that the generator can produce. */
+    static constexpr result_type (max)() noexcept
+    {
+        return (std::numeric_limits<result_type>::max)();
+    }
+
+    /** Returns the smallest value that the generator can produce. */
+    static constexpr result_type (min)() noexcept
+    {
+        return (std::numeric_limits<result_type>::min)();
+    }
+
+    /** Returns the state of the generator: the digits, oldest first, and then the carry. */
+    inline std::array<Digit, K> state() const noexcept
+    {
+        return state_;
+    }
+};
+
+} // namespace detail
+} // namespace random
+} // namespace boost
+
+#endif // BOOST_RANDOM_DETAIL_MWCX_BASE_HPP
diff --git a/include/boost/random/mwcx.hpp b/include/boost/random/mwcx.hpp
new file mode 100644
index 000000000..dcb0aac5a
--- /dev/null
+++ b/include/boost/random/mwcx.hpp
@@ -0,0 +1,288 @@
+/*
+ * Copyright Sebastiano Vigna 2026.
+ * Distributed under the Boost Software License, Version 1.0. (See
+ * accompanying file LICENSE_1_0.txt or copy at
+ * http://www.boost.org/LICENSE_1_0.txt)
+ *
+ * See http://www.boost.org for most recent version including documentation.
+ */
+
+#ifndef BOOST_RANDOM_MWCX_HPP
+#define BOOST_RANDOM_MWCX_HPP
+
+#include <boost/random/detail/config.hpp>
+#include <boost/random/detail/mwcx_arithmetic.hpp>
+#include <boost/random/detail/mwcx_base.hpp>
+#include <array>
+#include <cstdint>
+
+namespace boost {
+namespace random {
+
+/**
+ * This is MWC192X, a scrambled multiply-with-carry generator with 64-bit
+ * digits, lag two, and period approximately 2^191, from "Scrambled
+ * Multiply-with-Carry Pseudorandom Number Generators" by Sebastiano Vigna
+ * (https://arxiv.org/abs/2610.27313). It is recommended for all purposes:
+ * it is one of the fastest generators available and it passes all tests we
+ * are aware of. Its only 128-bit operations are a multiplication and a sum,
+ * and, since the lag is two, each multiplication uses a digit computed two
+ * steps before, so the processor can overlap the multiplications of
+ * consecutive steps.
+ *
+ * The output is scrambled by X: it is the most recent digit xored with the
+ * carry, as in David Thomas's MWC64X generator. The underlying generator
+ * simulates a multiplicative congruential generator with prime modulus
+ * m = 0xff4e24dd418f5433 * 2^128 - 1 and multiplier 2^(-64) mod m.
+ *
+ * discard() computes the state up to 2^64 - 1 steps ahead in time
+ * logarithmic in the number of steps, and jump(c, e) advances the generator
+ * by c * 2^e steps, which makes it possible to obtain non-overlapping
+ * subsequences at distances larger than discard() can express.
+ * jump(1, 96) is equivalent to the jump() function of the reference
+ * implementation: it can be used to generate 2^95 non-overlapping
+ * subsequences for parallel computations. jump(1, 144) is equivalent to
+ * long_jump(): it can be used to generate 2^47 starting points, from
+ * each of which jump(1, 96) will generate 2^48 non-overlapping
+ * subsequences for parallel distributed computations.
+ *
+ * The state is formed by two 64-bit digits x and y and a 64-bit carry c, and
+ * it is read as the integer x + y * 2^64 + c * 2^128, which must be nonzero
+ * and smaller than m. This generator is not suitable for cryptographic
+ * purposes.
+ */
+class mwc192x final : public detail::mwcx_base<mwc192x, std::uint64_t, 3, 3>
+{
+private:
+
+    using Base = detail::mwcx_base<mwc192x, std::uint64_t, 3, 3>;
+    friend Base;
+
+    static constexpr std::uint64_t a = UINT64_C(0xff4e24dd418f5433);
+
+    static words_type modulus() noexcept
+    {
+        return {{ UINT64_C(0xffffffffffffffff), UINT64_C(0xffffffffffffffff), a - 1 }};
+    }
+
+    static words_type inverse_base() noexcept
+    {
+        return {{ 0, a, 0 }};
+    }
+
+
+public:
+
+    using Base::Base;
+
+    inline result_type next() noexcept
+    {
+        const std::uint64_t result = state_[1] ^ state_[2];
+        std::uint64_t hi;
+        const std::uint64_t lo = detail::mwcx_mul_add(a, state_[0], state_[2], hi);
+        state_[0] = state_[1];
+        state_[1] = lo;
+        state_[2] = hi;
+
+        return result;
+    }
+};
+
+/**
+ * This is MWC128X, a scrambled multiply-with-carry generator with 64-bit
+ * digits, lag one, and period approximately 2^127, from "Scrambled
+ * Multiply-with-Carry Pseudorandom Number Generators" by Sebastiano Vigna
+ * (https://arxiv.org/abs/2610.27313). It is an alternative to mwc192x with
+ * a smaller state (128 bits), and it passes all tests we are aware of. It is
+ * slightly slower than mwc192x, as each multiplication needs the result of
+ * the previous step.
+ *
+ * The output is scrambled by X: it is the most recent digit xored with the
+ * carry, as in David Thomas's MWC64X generator. The underlying generator
+ * simulates a multiplicative congruential generator with prime modulus
+ * m = 0xff6da236274e252e * 2^64 - 1 and multiplier 2^(-64) mod m.
+ *
+ * discard() computes the state up to 2^64 - 1 steps ahead in time
+ * logarithmic in the number of steps, and jump(c, e) advances the generator
+ * by c * 2^e steps, which makes it possible to obtain non-overlapping
+ * subsequences at distances larger than discard() can express.
+ * jump(1, 64) is equivalent to the jump() function of the reference
+ * implementation: it can be used to generate 2^63 non-overlapping
+ * subsequences for parallel computations. jump(1, 96) is equivalent to
+ * long_jump(): it can be used to generate 2^31 starting points, from
+ * each of which jump(1, 64) will generate 2^32 non-overlapping
+ * subsequences for parallel distributed computations.
+ *
+ * The state is formed by a 64-bit digit x and a 64-bit carry c, and it is
+ * read as the integer x + c * 2^64, which must be nonzero and smaller than m.
+ * This generator is not suitable for cryptographic purposes.
+ */
+class mwc128x final : public detail::mwcx_base<mwc128x, std::uint64_t, 2, 2>
+{
+private:
+
+    using Base = detail::mwcx_base<mwc128x, std::uint64_t, 2, 2>;
+    friend Base;
+
+    static constexpr std::uint64_t a = UINT64_C(0xff6da236274e252e);
+
+    static words_type modulus() noexcept
+    {
+        return {{ UINT64_C(0xffffffffffffffff), a - 1 }};
+    }
+
+    static words_type inverse_base() noexcept
+    {
+        return {{ a, 0 }};
+    }
+
+
+public:
+
+    using Base::Base;
+
+    inline result_type next() noexcept
+    {
+        const std::uint64_t result = state_[0] ^ state_[1];
+        std::uint64_t hi;
+        state_[0] = detail::mwcx_mul_add(a, state_[0], state_[1], hi);
+        state_[1] = hi;
+
+        return result;
+    }
+};
+
+/**
+ * This is MWC128X/32, a scrambled multiply-with-carry generator with 32-bit
+ * digits, lag three, and period approximately 2^127, from "Scrambled
+ * Multiply-with-Carry Pseudorandom Number Generators" by Sebastiano Vigna
+ * (https://arxiv.org/abs/2610.27313). It returns 32-bit values, it is
+ * recommended for all purposes, and it passes all tests we are aware of. It
+ * needs just a 32 x 32 -> 64-bit multiplication.
+ *
+ * The output is scrambled by X: it is the most recent digit xored with the
+ * carry, as in David Thomas's MWC64X generator. The underlying generator
+ * simulates a multiplicative congruential generator with prime modulus
+ * m = 0xffcccc66 * 2^96 - 1 and multiplier 2^(-32) mod m.
+ *
+ * discard() computes the state up to 2^64 - 1 steps ahead in time
+ * logarithmic in the number of steps, and jump(c, e) advances the generator
+ * by c * 2^e steps, which makes it possible to obtain non-overlapping
+ * subsequences at distances larger than discard() can express.
+ * jump(1, 64) is equivalent to the jump() function of the reference
+ * implementation: it can be used to generate 2^63 non-overlapping
+ * subsequences for parallel computations. jump(1, 96) is equivalent to
+ * long_jump(): it can be used to generate 2^31 starting points, from
+ * each of which jump(1, 64) will generate 2^32 non-overlapping
+ * subsequences for parallel distributed computations.
+ *
+ * The state is formed by three 32-bit digits x, y, z and a 32-bit carry c,
+ * and it is read as the integer x + y * 2^32 + z * 2^64 + c * 2^96, which
+ * must be nonzero and smaller than m. This generator is not suitable for
+ * cryptographic purposes.
+ */
+class mwc128x32 final : public detail::mwcx_base<mwc128x32, std::uint32_t, 4, 2>
+{
+private:
+
+    using Base = detail::mwcx_base<mwc128x32, std::uint32_t, 4, 2>;
+    friend Base;
+
+    static constexpr std::uint32_t a = UINT32_C(0xffcccc66);
+
+    static words_type modulus() noexcept
+    {
+        return {{ UINT64_C(0xffffffffffffffff), static_cast<std::uint64_t>(a - 1) << 32 | UINT64_C(0xffffffff) }};
+    }
+
+    static words_type inverse_base() noexcept
+    {
+        return {{ 0, a }};
+    }
+
+
+public:
+
+    using Base::Base;
+
+    inline result_type next() noexcept
+    {
+        const std::uint32_t result = state_[2] ^ state_[3];
+        const std::uint64_t t = static_cast<std::uint64_t>(a) * state_[0] + state_[3];
+        state_[0] = state_[1];
+        state_[1] = state_[2];
+        state_[2] = static_cast<std::uint32_t>(t);
+        state_[3] = static_cast<std::uint32_t>(t >> 32);
+
+        return result;
+    }
+};
+
+/**
+ * This is MWC64X/32, a scrambled multiply-with-carry generator with 32-bit
+ * digits, lag one, and period approximately 2^63, from "Scrambled
+ * Multiply-with-Carry Pseudorandom Number Generators" by Sebastiano Vigna
+ * (https://arxiv.org/abs/2610.27313). It returns 32-bit values, and it is an
+ * alternative to mwc128x32 with a smaller state (64 bits), meant for
+ * environments in which registers are scarce, such as GPUs. It needs just a
+ * 32 x 32 -> 64-bit multiplication.
+ *
+ * The output is scrambled by X: it is the most recent digit xored with the
+ * carry, as in David Thomas's MWC64X generator. The underlying generator
+ * simulates a multiplicative congruential generator with prime modulus
+ * m = 0xfff145f8 * 2^32 - 1 and multiplier 2^(-32) mod m.
+ *
+ * discard() computes the state up to 2^64 - 1 steps ahead in time
+ * logarithmic in the number of steps, and jump(c, e) advances the generator
+ * by c * 2^e steps, which makes it possible to obtain non-overlapping
+ * subsequences at distances larger than discard() can express.
+ * jump(1, 32) is equivalent to the jump() function of the reference
+ * implementation: it can be used to generate 2^31 non-overlapping
+ * subsequences for parallel computations. jump(1, 48) is equivalent to
+ * long_jump(): it can be used to generate 2^15 starting points, from
+ * each of which jump(1, 32) will generate 2^16 non-overlapping
+ * subsequences for parallel distributed computations.
+ *
+ * The state is formed by a 32-bit digit x and a 32-bit carry c, and it is
+ * read as the integer x + c * 2^32, which must be nonzero and smaller than m.
+ * This generator is not suitable for cryptographic purposes.
+ */
+class mwc64x32 final : public detail::mwcx_base<mwc64x32, std::uint32_t, 2, 1>
+{
+private:
+
+    using Base = detail::mwcx_base<mwc64x32, std::uint32_t, 2, 1>;
+    friend Base;
+
+    static constexpr std::uint32_t a = UINT32_C(0xfff145f8);
+
+    static words_type modulus() noexcept
+    {
+        return {{ (static_cast<std::uint64_t>(a) << 32) - 1 }};
+    }
+
+    static words_type inverse_base() noexcept
+    {
+        return {{ a }};
+    }
+
+
+public:
+
+    using Base::Base;
+
+    inline result_type next() noexcept
+    {
+        const std::uint32_t result = state_[0] ^ state_[1];
+        const std::uint64_t t = static_cast<std::uint64_t>(a) * state_[0] + state_[1];
+        state_[0] = static_cast<std::uint32_t>(t);
+        state_[1] = static_cast<std::uint32_t>(t >> 32);
+
+        return result;
+    }
+};
+
+} // namespace random
+} // namespace boost
+
+#endif // BOOST_RANDOM_MWCX_HPP
diff --git a/performance/random_speed.cpp b/performance/random_speed.cpp
index 467e9a357..1b3df828e 100644
--- a/performance/random_speed.cpp
+++ b/performance/random_speed.cpp
@@ -392,6 +392,10 @@ int main(int argc, char*argv[])
   run(iter, "xoshiro128pp", boost::random::xoshiro128pp());
   run(iter, "xoshiro128f", boost::random::xoshiro128f());
   run(iter, "xoshiro128mm", boost::random::xoshiro128mm());
+  run(iter, "mwc192x", boost::random::mwc192x());
+  run(iter, "mwc128x", boost::random::mwc128x());
+  run(iter, "mwc128x32", boost::random::mwc128x32());
+  run(iter, "mwc64x32", boost::random::mwc64x32());
 
 #ifdef HAVE_MT19937INT_C
   // requires the original mt19937int.c
diff --git a/test/Jamfile.v2 b/test/Jamfile.v2
index 414b1c2a7..b8fcc5799 100644
--- a/test/Jamfile.v2
+++ b/test/Jamfile.v2
@@ -85,6 +85,18 @@ run test_xoshiro128mm.cpp /boost/test//boost_unit_test_framework ;
 run test_comp_xoshiro128mm.cpp ;
 run test_xoshiro128f.cpp /boost/test//boost_unit_test_framework ;
 run test_comp_xoshiro128f.cpp ;
+run test_mwc192x.cpp /boost/test//boost_unit_test_framework ;
+run test_comp_mwc192x.cpp ;
+run test_comp_mwc192x_portable.cpp ;
+run test_mwc128x.cpp /boost/test//boost_unit_test_framework ;
+run test_comp_mwc128x.cpp ;
+run test_comp_mwc128x_portable.cpp ;
+run test_mwc128x32.cpp /boost/test//boost_unit_test_framework ;
+run test_comp_mwc128x32.cpp ;
+run test_mwc64x32.cpp /boost/test//boost_unit_test_framework ;
+run test_comp_mwc64x32.cpp ;
+run test_mwcx_mul_add.cpp ;
+run test_mwcx_mul_add_portable.cpp ;
 
 run github_issue_133.cpp ;
 run github_issue_147.cpp ;
diff --git a/test/test_comp_mwc128x.cpp b/test/test_comp_mwc128x.cpp
new file mode 100644
index 000000000..493ff16c4
--- /dev/null
+++ b/test/test_comp_mwc128x.cpp
@@ -0,0 +1,161 @@
+/*
+ * Copyright Sebastiano Vigna 2026.
+ * Distributed under the Boost Software License, Version 1.0. (See
+ * accompanying file LICENSE_1_0.txt or copy at
+ * http://www.boost.org/LICENSE_1_0.txt)
+ *
+ * This file copies and pastes the original code for comparison under the following license
+ *
+ * Written in 2023-2026 by Sebastiano Vigna (vigna@acm.org)
+ *
+ * To the extent possible under law, the author has dedicated all copyright
+ * and related and neighboring rights to this software to the public domain
+ * worldwide.
+ *
+ * Permission to use, copy, modify, and/or distribute this software for any
+ * purpose with or without fee is hereby granted.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR
+ * IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+#include "test_comp_mwcx.hpp"
+#include <cstdint>
+
+using boost::random::mwc128x;
+using std::uint64_t;
+
+#if defined(__SIZEOF_INT128__)
+
+/* The reference implementation (https://github.com/vigna/MWCPRNG/blob/main/prngs/MWC128X.c),
+   with the state in a structure. */
+
+#define MWC_A1 0xff6da236274e252e
+
+struct reference
+{
+    uint64_t x, c;
+
+    uint64_t next(void) {
+        const uint64_t result = x ^ c;
+        const __uint128_t t = MWC_A1 * (__uint128_t)x + c;
+        x = t;
+        c = t >> 64;
+        return result;
+    }
+};
+
+static void compare_with_reference(const uint64_t x, const uint64_t c)
+{
+    reference r = { x, c };
+    mwc128x e = make_engine<mwc128x>({ x, c });
+    // The index of the first output that differs from the reference implementation
+    int i = 0;
+    while (i < 1000000 && e() == r.next())
+    {
+        ++i;
+    }
+
+    BOOST_TEST_EQ(i, 1000000);
+}
+
+#endif
+
+int main()
+{
+    const uint64_t a = UINT64_C(0xff6da236274e252e);
+    const uint64_t max = UINT64_C(0xffffffffffffffff);
+
+    // x = 0x0123456789abcdef, c = 1
+    const mwc128x start = make_engine<mwc128x>({ UINT64_C(0x0123456789abcdef), 1 });
+
+    // These values were produced with the reference implementation
+    check_outputs(start, { UINT64_C(0x0123456789abcdee), UINT64_C(0xba17a27d87b42a4a), UINT64_C(0xa622a4b120a9f538),
+                           UINT64_C(0x1bb7f645c6670559), UINT64_C(0x18529bad6097d72b), UINT64_C(0xf501042c00493549),
+                           UINT64_C(0x37e0328b3ce80699), UINT64_C(0x958cd94f7eb2e02f), UINT64_C(0x68e32730e418e30a),
+                           UINT64_C(0x4519e7cc8909e77d) });
+
+    // discard() at large distances, with states computed by the reference implementation
+    mwc128x e;
+    e = start;
+    e.discard(12345);
+    BOOST_TEST(e == make_engine<mwc128x>({ UINT64_C(0xe7ef89c7c08b0a84), UINT64_C(0x63f3b10034fa793a) }));
+    e = start;
+    e.discard(UINT64_C(0x9e3779b97f4a7c15));
+    BOOST_TEST(e == make_engine<mwc128x>({ UINT64_C(0x7f0a05c03102a9ac), UINT64_C(0xe4cf2f5a527eb8a2) }));
+    e = start;
+    e.discard(UINT64_C(0xffffffffffffffff));
+    BOOST_TEST(e == make_engine<mwc128x>({ UINT64_C(0xe4f3ec1d416d5dad), UINT64_C(0x57bae73c39fc0c2f) }));
+    e = start;
+    e.discard(UINT64_C(1) << 32);
+    BOOST_TEST(e == make_engine<mwc128x>({ UINT64_C(0xbdfa9e1fd955e762), UINT64_C(0x1347546687e55513) }));
+    e = start;
+    e.discard(UINT64_C(1) << 48);
+    BOOST_TEST(e == make_engine<mwc128x>({ UINT64_C(0x3b18b6c941ea29e2), UINT64_C(0x8eacc536be1701a0) }));
+    // 2^64, the distance of jump() in the reference implementation
+    e = start;
+    e.discard(UINT64_C(1) << 63);
+    e.discard(UINT64_C(1) << 63);
+    BOOST_TEST(e == make_engine<mwc128x>({ UINT64_C(0x5a936408afe2e245), UINT64_C(0xe471051f76a629fe) }));
+
+    // jump(c, e), with states computed by the reference implementation:
+    // jump(1, 64) and jump(1, 96) are its jump() and long_jump()
+    e = start;
+    e.jump(1, 64);
+    BOOST_TEST(e == make_engine<mwc128x>({ UINT64_C(0x5a936408afe2e245), UINT64_C(0xe471051f76a629fe) }));
+    e = start;
+    e.jump(1, 96);
+    BOOST_TEST(e == make_engine<mwc128x>({ UINT64_C(0x7e2cba541c5f4c51), UINT64_C(0x30a344ed8d373892) }));
+    e = start;
+    e.jump(3, 70);
+    BOOST_TEST(e == make_engine<mwc128x>({ UINT64_C(0x34d3926d1a7906f8), UINT64_C(0x131cbd0e516ab608) }));
+    e = start;
+    e.jump(12345, 0);
+    BOOST_TEST(e == make_engine<mwc128x>({ UINT64_C(0xe7ef89c7c08b0a84), UINT64_C(0x63f3b10034fa793a) }));
+
+    // discard() from the default state, a distance computed with the reference implementation
+    e = mwc128x();
+    e.discard(1000000);
+    BOOST_TEST(e == make_engine<mwc128x>({ UINT64_C(0x45e7e0993ad51097), UINT64_C(0xaad79e0ed8884d55) }));
+
+    check_discard(start);
+    check_discard(mwc128x(42));
+    check_discard_composition(start);
+    check_jump(start, 64);
+
+    check_streaming(start);
+    check_short_seed<mwc128x>();
+
+    // Seeding with an integer: SplitMix digit and carry one
+    {
+        uint64_t s = 42;
+        uint64_t w = (s += UINT64_C(0x9E3779B97F4A7C15));
+        w = (w ^ (w >> 30)) * UINT64_C(0xBF58476D1CE4E5B9);
+        w = (w ^ (w >> 27)) * UINT64_C(0x94D049BB133111EB);
+        w ^= w >> 31;
+        BOOST_TEST(mwc128x(42) == make_engine<mwc128x>({ w, 1 }));
+        BOOST_TEST(mwc128x() == mwc128x(0));
+    }
+
+    // The zero state and the states not smaller than m get carry one; the
+    // other states, including the largest one, are left unchanged.
+    BOOST_TEST(make_engine<mwc128x>({ 0, 0 }).state()[1] == 1);
+    BOOST_TEST(make_engine<mwc128x>({ max, a - 1 }).state()[1] == 1);
+    BOOST_TEST(make_engine<mwc128x>({ 5, max }).state()[1] == 1);
+    BOOST_TEST(make_engine<mwc128x>({ max - 1, a - 1 }).state()[1] == a - 1);
+    BOOST_TEST(make_engine<mwc128x>({ 1, 0 }).state()[1] == 0);
+
+#if defined(__SIZEOF_INT128__)
+    compare_with_reference(UINT64_C(0x0123456789abcdef), 1);
+    compare_with_reference(max - 1, a - 1);
+    compare_with_reference(max, a - 2);
+    compare_with_reference(1, 0);
+#endif
+
+    return boost::report_errors();
+}
diff --git a/test/test_comp_mwc128x32.cpp b/test/test_comp_mwc128x32.cpp
new file mode 100644
index 000000000..42217dd66
--- /dev/null
+++ b/test/test_comp_mwc128x32.cpp
@@ -0,0 +1,162 @@
+/*
+ * Copyright Sebastiano Vigna 2026.
+ * Distributed under the Boost Software License, Version 1.0. (See
+ * accompanying file LICENSE_1_0.txt or copy at
+ * http://www.boost.org/LICENSE_1_0.txt)
+ *
+ * This file copies and pastes the original code for comparison under the following license
+ *
+ * Written in 2023-2026 by Sebastiano Vigna (vigna@acm.org)
+ *
+ * To the extent possible under law, the author has dedicated all copyright
+ * and related and neighboring rights to this software to the public domain
+ * worldwide.
+ *
+ * Permission to use, copy, modify, and/or distribute this software for any
+ * purpose with or without fee is hereby granted.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR
+ * IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+#include "test_comp_mwcx.hpp"
+#include <cstdint>
+
+using boost::random::mwc128x32;
+using std::uint32_t;
+using std::uint64_t;
+
+/* The reference implementation (https://github.com/vigna/MWCPRNG/blob/main/prngs/MWC128X_32.c),
+   with the state in a structure. */
+
+#define MWC_A3 0xffcccc66
+
+struct reference
+{
+    uint32_t x, y, z, c;
+
+    uint32_t next(void) {
+        const uint32_t result = z ^ c;
+        const uint64_t t = MWC_A3 * (uint64_t)x + c;
+        x = y;
+        y = z;
+        z = t;
+        c = t >> 32;
+        return result;
+    }
+};
+
+static void compare_with_reference(const uint32_t x, const uint32_t y, const uint32_t z, const uint32_t c)
+{
+    reference r = { x, y, z, c };
+    mwc128x32 e = make_engine<mwc128x32>({ x, y, z, c });
+    // The index of the first output that differs from the reference implementation
+    int i = 0;
+    while (i < 1000000 && e() == r.next())
+    {
+        ++i;
+    }
+
+    BOOST_TEST_EQ(i, 1000000);
+}
+
+int main()
+{
+    const uint32_t a = UINT32_C(0xffcccc66);
+    const uint32_t max = UINT32_C(0xffffffff);
+
+    // x = 0x89abcdef, y = 0x76543210, z = 0x4b5a6978, c = 1
+    const mwc128x32 start = make_engine<mwc128x32>({ UINT32_C(0x89abcdef), UINT32_C(0x76543210), UINT32_C(0x4b5a6978), 1 });
+
+    // These values were produced with the reference implementation
+    check_outputs(start, { UINT32_C(0x4b5a6979), UINT32_C(0x6792c5cf), UINT32_C(0xedfc7024), UINT32_C(0x50a57a05),
+                           UINT32_C(0x329b3e06), UINT32_C(0x1f1db247), UINT32_C(0x4c7e40d2), UINT32_C(0x3cbaba78),
+                           UINT32_C(0x687c6b34), UINT32_C(0xe8c71e6d) });
+
+    // discard() at large distances, with states computed by the reference implementation
+    mwc128x32 e;
+    e = start;
+    e.discard(12345);
+    BOOST_TEST(e == make_engine<mwc128x32>({ UINT32_C(0x2602fd2c), UINT32_C(0xd27ae5a4), UINT32_C(0x6dc5ec3d), UINT32_C(0x593ac844) }));
+    e = start;
+    e.discard(UINT64_C(0x9e3779b97f4a7c15));
+    BOOST_TEST(e == make_engine<mwc128x32>({ UINT32_C(0x6385879f), UINT32_C(0xf62de3b4), UINT32_C(0x3028e545), UINT32_C(0xf801572e) }));
+    e = start;
+    e.discard(UINT64_C(0xffffffffffffffff));
+    BOOST_TEST(e == make_engine<mwc128x32>({ UINT32_C(0x4e12e1e3), UINT32_C(0x69ec7c23), UINT32_C(0x85620075), UINT32_C(0x6071af3a) }));
+    e = start;
+    e.discard(UINT64_C(1) << 32);
+    BOOST_TEST(e == make_engine<mwc128x32>({ UINT32_C(0x3988a8aa), UINT32_C(0xe38a592d), UINT32_C(0x18281c8a), UINT32_C(0xfd2a52db) }));
+    e = start;
+    e.discard(UINT64_C(1) << 48);
+    BOOST_TEST(e == make_engine<mwc128x32>({ UINT32_C(0x2b1169fc), UINT32_C(0x28df4959), UINT32_C(0x08847838), UINT32_C(0x3d0297fd) }));
+    // 2^64, the distance of jump() in the reference implementation
+    e = start;
+    e.discard(UINT64_C(1) << 63);
+    e.discard(UINT64_C(1) << 63);
+    BOOST_TEST(e == make_engine<mwc128x32>({ UINT32_C(0x69ec7c23), UINT32_C(0x85620075), UINT32_C(0xa5dc93ac), UINT32_C(0x4e034463) }));
+
+    // jump(c, e), with states computed by the reference implementation:
+    // jump(1, 64) and jump(1, 96) are its jump() and long_jump()
+    e = start;
+    e.jump(1, 64);
+    BOOST_TEST(e == make_engine<mwc128x32>({ UINT32_C(0x69ec7c23), UINT32_C(0x85620075), UINT32_C(0xa5dc93ac), UINT32_C(0x4e034463) }));
+    e = start;
+    e.jump(1, 96);
+    BOOST_TEST(e == make_engine<mwc128x32>({ UINT32_C(0xee6eebbd), UINT32_C(0x140f21ca), UINT32_C(0x45518300), UINT32_C(0x1321a37d) }));
+    e = start;
+    e.jump(3, 70);
+    BOOST_TEST(e == make_engine<mwc128x32>({ UINT32_C(0x8cb9387d), UINT32_C(0x68a3bfad), UINT32_C(0x2d7e249b), UINT32_C(0x7b14417a) }));
+    e = start;
+    e.jump(12345, 0);
+    BOOST_TEST(e == make_engine<mwc128x32>({ UINT32_C(0x2602fd2c), UINT32_C(0xd27ae5a4), UINT32_C(0x6dc5ec3d), UINT32_C(0x593ac844) }));
+
+    // discard() from the default state, a distance computed with the reference implementation
+    e = mwc128x32();
+    e.discard(1000000);
+    BOOST_TEST(e == make_engine<mwc128x32>({ UINT32_C(0x67dd14bf), UINT32_C(0x2079e61c), UINT32_C(0x08abc63e), UINT32_C(0xcbd6c060) }));
+
+    check_discard(start);
+    check_discard(mwc128x32(42));
+    check_discard_composition(start);
+    check_jump(start, 64);
+
+    check_streaming(start);
+    check_short_seed<mwc128x32>();
+
+    // Seeding with an integer: the lower and upper halves of two SplitMix
+    // outputs (the upper half of the second is not used), and carry one
+    {
+        uint64_t s = 42;
+        uint64_t z[2];
+        for (auto& w : z)
+        {
+            w = (s += UINT64_C(0x9E3779B97F4A7C15));
+            w = (w ^ (w >> 30)) * UINT64_C(0xBF58476D1CE4E5B9);
+            w = (w ^ (w >> 27)) * UINT64_C(0x94D049BB133111EB);
+            w ^= w >> 31;
+        }
+        BOOST_TEST(mwc128x32(42) == make_engine<mwc128x32>({ static_cast<uint32_t>(z[0]), static_cast<uint32_t>(z[0] >> 32), static_cast<uint32_t>(z[1]), 1 }));
+        BOOST_TEST(mwc128x32() == mwc128x32(0));
+    }
+
+    // The zero state and the states not smaller than m get carry one; the
+    // other states, including the largest one, are left unchanged.
+    BOOST_TEST(make_engine<mwc128x32>({ 0, 0, 0, 0 }).state()[3] == 1);
+    BOOST_TEST(make_engine<mwc128x32>({ max, max, max, a - 1 }).state()[3] == 1);
+    BOOST_TEST(make_engine<mwc128x32>({ 5, 0, 0, max }).state()[3] == 1);
+    BOOST_TEST(make_engine<mwc128x32>({ max, max, max - 1, a - 1 }).state()[3] == a - 1);
+    BOOST_TEST(make_engine<mwc128x32>({ 1, 0, 0, 0 }).state()[3] == 0);
+
+    compare_with_reference(UINT32_C(0x89abcdef), UINT32_C(0x76543210), UINT32_C(0x4b5a6978), 1);
+    compare_with_reference(max, max, max - 1, a - 1);
+    compare_with_reference(max, max, max, a - 2);
+    compare_with_reference(1, 0, 0, 0);
+
+    return boost::report_errors();
+}
diff --git a/test/test_comp_mwc128x_portable.cpp b/test/test_comp_mwc128x_portable.cpp
new file mode 100644
index 000000000..4051d3f9b
--- /dev/null
+++ b/test/test_comp_mwc128x_portable.cpp
@@ -0,0 +1,11 @@
+/*
+ * Copyright Sebastiano Vigna 2026.
+ * Distributed under the Boost Software License, Version 1.0. (See
+ * accompanying file LICENSE_1_0.txt or copy at
+ * http://www.boost.org/LICENSE_1_0.txt)
+ *
+ * The tests of test_comp_mwc128x.cpp with the portable 64 x 64 -> 128-bit multiplication.
+ */
+
+#define BOOST_RANDOM_MWCX_PORTABLE_MUL
+#include "test_comp_mwc128x.cpp"
diff --git a/test/test_comp_mwc192x.cpp b/test/test_comp_mwc192x.cpp
new file mode 100644
index 000000000..d54c719f2
--- /dev/null
+++ b/test/test_comp_mwc192x.cpp
@@ -0,0 +1,166 @@
+/*
+ * Copyright Sebastiano Vigna 2026.
+ * Distributed under the Boost Software License, Version 1.0. (See
+ * accompanying file LICENSE_1_0.txt or copy at
+ * http://www.boost.org/LICENSE_1_0.txt)
+ *
+ * This file copies and pastes the original code for comparison under the following license
+ *
+ * Written in 2023-2026 by Sebastiano Vigna (vigna@acm.org)
+ *
+ * To the extent possible under law, the author has dedicated all copyright
+ * and related and neighboring rights to this software to the public domain
+ * worldwide.
+ *
+ * Permission to use, copy, modify, and/or distribute this software for any
+ * purpose with or without fee is hereby granted.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR
+ * IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+#include "test_comp_mwcx.hpp"
+#include <cstdint>
+
+using boost::random::mwc192x;
+using std::uint64_t;
+
+#if defined(__SIZEOF_INT128__)
+
+/* The reference implementation (https://github.com/vigna/MWCPRNG/blob/main/prngs/MWC192X.c),
+   with the state in a structure. */
+
+#define MWC_A2 0xff4e24dd418f5433
+
+struct reference
+{
+    uint64_t x, y, c;
+
+    uint64_t next(void) {
+        const uint64_t result = y ^ c;
+        const __uint128_t t = MWC_A2 * (__uint128_t)x + c;
+        x = y;
+        y = t;
+        c = t >> 64;
+        return result;
+    }
+};
+
+static void compare_with_reference(const uint64_t x, const uint64_t y, const uint64_t c)
+{
+    reference r = { x, y, c };
+    mwc192x e = make_engine<mwc192x>({ x, y, c });
+    // The index of the first output that differs from the reference implementation
+    int i = 0;
+    while (i < 1000000 && e() == r.next())
+    {
+        ++i;
+    }
+
+    BOOST_TEST_EQ(i, 1000000);
+}
+
+#endif
+
+int main()
+{
+    const uint64_t a = UINT64_C(0xff4e24dd418f5433);
+    const uint64_t max = UINT64_C(0xffffffffffffffff);
+
+    // x = 0x0123456789abcdef, y = 0xfedcba9876543210, c = 1
+    const mwc192x start = make_engine<mwc192x>({ UINT64_C(0x0123456789abcdef), UINT64_C(0xfedcba9876543210), 1 });
+
+    // These values were produced with the reference implementation
+    check_outputs(start, { UINT64_C(0xfedcba9876543211), UINT64_C(0xac8dcd6abd635774), UINT64_C(0xaa0f361e7b30715d),
+                           UINT64_C(0xe104b15c76fca2e4), UINT64_C(0x2c12e1a7db8cbfa1), UINT64_C(0x6105bf55a327d362),
+                           UINT64_C(0x62ace12dcd4dc6ba), UINT64_C(0x90e3539226e6de34), UINT64_C(0x112e01165d99bb15),
+                           UINT64_C(0x607f209a25bea6ba) });
+
+    // discard() at large distances, with states computed by the reference implementation
+    mwc192x e;
+    e = start;
+    e.discard(12345);
+    BOOST_TEST(e == make_engine<mwc192x>({ UINT64_C(0x10ef1947237ccbb6), UINT64_C(0x4de156bc5cf48087), UINT64_C(0xe8c5b59ed41706bd) }));
+    e = start;
+    e.discard(UINT64_C(0x9e3779b97f4a7c15));
+    BOOST_TEST(e == make_engine<mwc192x>({ UINT64_C(0x9a9addfe20bc00c8), UINT64_C(0xbf7da39fa67b2851), UINT64_C(0x1adbc6886c40535a) }));
+    e = start;
+    e.discard(UINT64_C(0xffffffffffffffff));
+    BOOST_TEST(e == make_engine<mwc192x>({ UINT64_C(0x152d4c0143c9cfd7), UINT64_C(0x48eb25ea11450b01), UINT64_C(0x6b0f2e23d7b82181) }));
+    e = start;
+    e.discard(UINT64_C(1) << 32);
+    BOOST_TEST(e == make_engine<mwc192x>({ UINT64_C(0x93e9a4a82fd6439b), UINT64_C(0x1da0e1033c9c2767), UINT64_C(0x3c6f78083948680f) }));
+    e = start;
+    e.discard(UINT64_C(1) << 48);
+    BOOST_TEST(e == make_engine<mwc192x>({ UINT64_C(0x33df4ff1df935014), UINT64_C(0xb479d139117b3c1b), UINT64_C(0x042ebef1828ce592) }));
+    // 2^64
+    e = start;
+    e.discard(UINT64_C(1) << 63);
+    e.discard(UINT64_C(1) << 63);
+    BOOST_TEST(e == make_engine<mwc192x>({ UINT64_C(0x48eb25ea11450b01), UINT64_C(0xa9e580be41381556), UINT64_C(0x151e958f171fd3b5) }));
+
+    // jump(c, e), with states computed by the reference implementation:
+    // jump(1, 96) and jump(1, 144) are its jump() and long_jump()
+    e = start;
+    e.jump(1, 96);
+    BOOST_TEST(e == make_engine<mwc192x>({ UINT64_C(0x5438b996afe69147), UINT64_C(0x4d8061a4a1f2862f), UINT64_C(0xd600e9332b37e5a7) }));
+    e = start;
+    e.jump(1, 144);
+    BOOST_TEST(e == make_engine<mwc192x>({ UINT64_C(0x9ff4a43e2bc6a812), UINT64_C(0xf695800ae27947cf), UINT64_C(0x008f2e60cfea1794) }));
+    e = start;
+    e.jump(3, 70);
+    BOOST_TEST(e == make_engine<mwc192x>({ UINT64_C(0x746e2db78364cb25), UINT64_C(0x5b3243b447ce93d9), UINT64_C(0x6e64eea0ebc89f14) }));
+    e = start;
+    e.jump(12345, 0);
+    BOOST_TEST(e == make_engine<mwc192x>({ UINT64_C(0x10ef1947237ccbb6), UINT64_C(0x4de156bc5cf48087), UINT64_C(0xe8c5b59ed41706bd) }));
+
+    // discard() from the default state, a distance computed with the reference implementation
+    e = mwc192x();
+    e.discard(1000000);
+    BOOST_TEST(e == make_engine<mwc192x>({ UINT64_C(0xcd2e0328326ef24d), UINT64_C(0x01cf511a1e056113), UINT64_C(0xe263fb805f0f7665) }));
+
+    check_discard(start);
+    check_discard(mwc192x(42));
+    check_discard_composition(start);
+    check_jump(start, 96);
+
+    check_streaming(start);
+    check_short_seed<mwc192x>();
+
+    // Seeding with an integer: SplitMix digits and carry one
+    {
+        uint64_t s = 42;
+        uint64_t z[2];
+        for (auto& w : z)
+        {
+            w = (s += UINT64_C(0x9E3779B97F4A7C15));
+            w = (w ^ (w >> 30)) * UINT64_C(0xBF58476D1CE4E5B9);
+            w = (w ^ (w >> 27)) * UINT64_C(0x94D049BB133111EB);
+            w ^= w >> 31;
+        }
+        BOOST_TEST(mwc192x(42) == make_engine<mwc192x>({ z[0], z[1], 1 }));
+        BOOST_TEST(mwc192x() == mwc192x(0));
+    }
+
+    // The zero state and the states not smaller than m get carry one; the
+    // other states, including the largest one, are left unchanged.
+    BOOST_TEST(make_engine<mwc192x>({ 0, 0, 0 }).state()[2] == 1);
+    BOOST_TEST(make_engine<mwc192x>({ max, max, a - 1 }).state()[2] == 1);
+    BOOST_TEST(make_engine<mwc192x>({ 5, 0, max }).state()[2] == 1);
+    BOOST_TEST(make_engine<mwc192x>({ max, max - 1, a - 1 }).state()[2] == a - 1);
+    BOOST_TEST(make_engine<mwc192x>({ 1, 0, 0 }).state()[2] == 0);
+
+#if defined(__SIZEOF_INT128__)
+    compare_with_reference(UINT64_C(0x0123456789abcdef), UINT64_C(0xfedcba9876543210), 1);
+    compare_with_reference(max, max - 1, a - 1);
+    compare_with_reference(max, max, a - 2);
+    compare_with_reference(1, 0, 0);
+#endif
+
+    return boost::report_errors();
+}
diff --git a/test/test_comp_mwc192x_portable.cpp b/test/test_comp_mwc192x_portable.cpp
new file mode 100644
index 000000000..1f14b4f63
--- /dev/null
+++ b/test/test_comp_mwc192x_portable.cpp
@@ -0,0 +1,11 @@
+/*
+ * Copyright Sebastiano Vigna 2026.
+ * Distributed under the Boost Software License, Version 1.0. (See
+ * accompanying file LICENSE_1_0.txt or copy at
+ * http://www.boost.org/LICENSE_1_0.txt)
+ *
+ * The tests of test_comp_mwc192x.cpp with the portable 64 x 64 -> 128-bit multiplication.
+ */
+
+#define BOOST_RANDOM_MWCX_PORTABLE_MUL
+#include "test_comp_mwc192x.cpp"
diff --git a/test/test_comp_mwc64x32.cpp b/test/test_comp_mwc64x32.cpp
new file mode 100644
index 000000000..d8cc05e27
--- /dev/null
+++ b/test/test_comp_mwc64x32.cpp
@@ -0,0 +1,155 @@
+/*
+ * Copyright Sebastiano Vigna 2026.
+ * Distributed under the Boost Software License, Version 1.0. (See
+ * accompanying file LICENSE_1_0.txt or copy at
+ * http://www.boost.org/LICENSE_1_0.txt)
+ *
+ * This file copies and pastes the original code for comparison under the following license
+ *
+ * Written in 2023-2026 by Sebastiano Vigna (vigna@acm.org)
+ *
+ * To the extent possible under law, the author has dedicated all copyright
+ * and related and neighboring rights to this software to the public domain
+ * worldwide.
+ *
+ * Permission to use, copy, modify, and/or distribute this software for any
+ * purpose with or without fee is hereby granted.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR
+ * IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+#include "test_comp_mwcx.hpp"
+#include <cstdint>
+
+using boost::random::mwc64x32;
+using std::uint32_t;
+using std::uint64_t;
+
+/* The reference implementation (https://github.com/vigna/MWCPRNG/blob/main/prngs/MWC64X_32.c),
+   with the state in a structure. */
+
+#define MWC_A1 0xfff145f8
+
+struct reference
+{
+    uint32_t x, c;
+
+    uint32_t next(void) {
+        const uint32_t result = x ^ c;
+        const uint64_t t = MWC_A1 * (uint64_t)x + c;
+        x = t;
+        c = t >> 32;
+        return result;
+    }
+};
+
+static void compare_with_reference(const uint32_t x, const uint32_t c)
+{
+    reference r = { x, c };
+    mwc64x32 e = make_engine<mwc64x32>({ x, c });
+    // The index of the first output that differs from the reference implementation
+    int i = 0;
+    while (i < 1000000 && e() == r.next())
+    {
+        ++i;
+    }
+
+    BOOST_TEST_EQ(i, 1000000);
+}
+
+int main()
+{
+    const uint32_t a = UINT32_C(0xfff145f8);
+    const uint32_t max = UINT32_C(0xffffffff);
+
+    // x = 0x89abcdef, c = 1
+    const mwc64x32 start = make_engine<mwc64x32>({ UINT32_C(0x89abcdef), 1 });
+
+    // These values were produced with the reference implementation
+    check_outputs(start, { UINT32_C(0x89abcdee), UINT32_C(0x124c08f7), UINT32_C(0x1dd9f62f), UINT32_C(0xe3dcdf5f),
+                           UINT32_C(0x6aadea99), UINT32_C(0x3b18a0d9), UINT32_C(0x8a555e1f), UINT32_C(0x46d69f8e),
+                           UINT32_C(0x9edf4569), UINT32_C(0x85ccb634) });
+
+    // discard() at large distances, with states computed by the reference implementation
+    mwc64x32 e;
+    e = start;
+    e.discard(12345);
+    BOOST_TEST(e == make_engine<mwc64x32>({ UINT32_C(0x8304459c), UINT32_C(0xc3636383) }));
+    e = start;
+    e.discard(UINT64_C(0x9e3779b97f4a7c15));
+    BOOST_TEST(e == make_engine<mwc64x32>({ UINT32_C(0x56b9861b), UINT32_C(0x0a7a33b3) }));
+    e = start;
+    e.discard(UINT64_C(0xffffffffffffffff));
+    BOOST_TEST(e == make_engine<mwc64x32>({ UINT32_C(0x0b5d9549), UINT32_C(0x367c459d) }));
+    e = start;
+    e.discard(UINT64_C(1) << 32); // jump() of the reference implementation
+    BOOST_TEST(e == make_engine<mwc64x32>({ UINT32_C(0xf9132140), UINT32_C(0xd28b3488) }));
+    e = start;
+    e.discard(UINT64_C(1) << 48); // long_jump() of the reference implementation
+    BOOST_TEST(e == make_engine<mwc64x32>({ UINT32_C(0x48cd3df2), UINT32_C(0x5b545e15) }));
+    // 2^64
+    e = start;
+    e.discard(UINT64_C(1) << 63);
+    e.discard(UINT64_C(1) << 63);
+    BOOST_TEST(e == make_engine<mwc64x32>({ UINT32_C(0xb31a9155), UINT32_C(0x0b5cede8) }));
+
+    // jump(c, e), with states computed by the reference implementation:
+    // jump(1, 32) and jump(1, 48) are its jump() and long_jump()
+    e = start;
+    e.jump(1, 32);
+    BOOST_TEST(e == make_engine<mwc64x32>({ UINT32_C(0xf9132140), UINT32_C(0xd28b3488) }));
+    e = start;
+    e.jump(1, 48);
+    BOOST_TEST(e == make_engine<mwc64x32>({ UINT32_C(0x48cd3df2), UINT32_C(0x5b545e15) }));
+    e = start;
+    e.jump(3, 70);
+    BOOST_TEST(e == make_engine<mwc64x32>({ UINT32_C(0xf99803f0), UINT32_C(0x8856dbf7) }));
+    e = start;
+    e.jump(12345, 0);
+    BOOST_TEST(e == make_engine<mwc64x32>({ UINT32_C(0x8304459c), UINT32_C(0xc3636383) }));
+
+    // discard() from the default state, a distance computed with the reference implementation
+    e = mwc64x32();
+    e.discard(1000000);
+    BOOST_TEST(e == make_engine<mwc64x32>({ UINT32_C(0x0afc4416), UINT32_C(0x2a3869bb) }));
+
+    check_discard(start);
+    check_discard(mwc64x32(42));
+    check_discard_composition(start);
+    check_jump(start, 32);
+
+    check_streaming(start);
+    check_short_seed<mwc64x32>();
+
+    // Seeding with an integer: the lower half of a SplitMix output, and carry one
+    {
+        uint64_t s = 42;
+        uint64_t w = (s += UINT64_C(0x9E3779B97F4A7C15));
+        w = (w ^ (w >> 30)) * UINT64_C(0xBF58476D1CE4E5B9);
+        w = (w ^ (w >> 27)) * UINT64_C(0x94D049BB133111EB);
+        w ^= w >> 31;
+        BOOST_TEST(mwc64x32(42) == make_engine<mwc64x32>({ static_cast<uint32_t>(w), 1 }));
+        BOOST_TEST(mwc64x32() == mwc64x32(0));
+    }
+
+    // The zero state and the states not smaller than m get carry one; the
+    // other states, including the largest one, are left unchanged.
+    BOOST_TEST(make_engine<mwc64x32>({ 0, 0 }).state()[1] == 1);
+    BOOST_TEST(make_engine<mwc64x32>({ max, a - 1 }).state()[1] == 1);
+    BOOST_TEST(make_engine<mwc64x32>({ 5, max }).state()[1] == 1);
+    BOOST_TEST(make_engine<mwc64x32>({ max - 1, a - 1 }).state()[1] == a - 1);
+    BOOST_TEST(make_engine<mwc64x32>({ 1, 0 }).state()[1] == 0);
+
+    compare_with_reference(UINT32_C(0x89abcdef), 1);
+    compare_with_reference(max - 1, a - 1);
+    compare_with_reference(max, a - 2);
+    compare_with_reference(1, 0);
+
+    return boost::report_errors();
+}
diff --git a/test/test_comp_mwcx.hpp b/test/test_comp_mwcx.hpp
new file mode 100644
index 000000000..430d68e2f
--- /dev/null
+++ b/test/test_comp_mwcx.hpp
@@ -0,0 +1,141 @@
+/*
+ * Copyright Sebastiano Vigna 2026.
+ * Distributed under the Boost Software License, Version 1.0. (See
+ * accompanying file LICENSE_1_0.txt or copy at
+ * http://www.boost.org/LICENSE_1_0.txt)
+ *
+ * Common checks of the scrambled multiply-with-carry generators (mwcx.hpp).
+ */
+
+#ifndef BOOST_RANDOM_TEST_COMP_MWCX_HPP
+#define BOOST_RANDOM_TEST_COMP_MWCX_HPP
+
+#include <boost/random/mwcx.hpp>
+#include <boost/core/lightweight_test.hpp>
+#include <array>
+#include <cstdint>
+#include <initializer_list>
+#include <sstream>
+#include <stdexcept>
+#include <vector>
+
+// Builds a generator from its state: the digits, oldest first, and then the carry.
+template <typename Engine>
+Engine make_engine(std::initializer_list<typename Engine::result_type> state)
+{
+    const std::vector<typename Engine::result_type> v(state);
+    Engine e;
+    e.seed(v.begin(), v.end());
+    return e;
+}
+
+// Checks that the first outputs of a generator are the expected ones.
+template <typename Engine>
+void check_outputs(Engine e, std::initializer_list<typename Engine::result_type> expected)
+{
+    for (const auto x : expected)
+    {
+        BOOST_TEST_EQ(e(), x);
+    }
+}
+
+// Checks that discard(d) is equivalent to d calls to operator()() and to
+// jump(d, 0), across the thresholds below which discard() steps the generator.
+template <typename Engine>
+void check_discard(const Engine& start)
+{
+    static const std::uint64_t distances[] = { 0, 1, 2, 3, 7, 64, 1000, 4095, 4096, 4097, 9307, 16383, 16384, 16385, 65535, 65536, 65537, 100000 };
+    for (const std::uint64_t d : distances)
+    {
+        Engine a = start;
+        for (std::uint64_t i = 0; i < d; ++i)
+        {
+            a();
+        }
+
+        Engine b = start;
+        b.discard(d);
+        BOOST_TEST(a == b);
+
+        Engine c = start;
+        c.jump(d, 0);
+        BOOST_TEST(a == c);
+    }
+}
+
+// Checks that discard(u) followed by discard(v) is equivalent to discard(u + v).
+template <typename Engine>
+void check_discard_composition(const Engine& start)
+{
+    static const std::uint64_t pairs[][2] = {
+        { 1, UINT64_C(0xfffffffffffffffe) },
+        { 100000, UINT64_C(0x123456789abcdef) },
+        { UINT64_C(0x7fffffffffffffff), UINT64_C(0x8000000000000000) },
+        { UINT64_C(0x9e3779b97f4a7c15), UINT64_C(0x0123456789abcdef) }
+    };
+    for (const auto& p : pairs)
+    {
+        Engine a = start;
+        a.discard(p[0]);
+        a.discard(p[1]);
+        Engine b = start;
+        b.discard(p[0] + p[1]);
+        BOOST_TEST(a == b);
+    }
+}
+
+// Checks that jump(c, e) is equivalent to jump(c * 2^e, 0) when c * 2^e
+// fits in 64 bits, and that jump(c1, e) followed by jump(c2, e) is
+// equivalent to jump(c1 + c2, e).
+template <typename Engine>
+void check_jump(const Engine& start, const unsigned int jump_log)
+{
+    static const std::uint64_t ce[][2] = { { 3, 8 }, { 5, 20 }, { UINT64_C(0x1234567), 30 }, { 1, 63 }, { UINT64_C(0x9e3779b97f4a7c15), 0 } };
+    for (const auto& p : ce)
+    {
+        Engine a = start;
+        a.jump(p[0], static_cast<unsigned int>(p[1]));
+        Engine b = start;
+        b.jump(p[0] << p[1], 0);
+        BOOST_TEST(a == b);
+    }
+
+    const unsigned int es[] = { 0, 1, 64, jump_log, 2 * jump_log };
+    for (const unsigned int e : es)
+    {
+        Engine a = start;
+        a.jump(3, e);
+        a.jump(UINT64_C(0x123456789abcdef), e);
+        Engine b = start;
+        b.jump(3 + UINT64_C(0x123456789abcdef), e);
+        BOOST_TEST(a == b);
+    }
+}
+
+// Checks that operator<< and operator>> preserve the state.
+template <typename Engine>
+void check_streaming(const Engine& e)
+{
+    std::stringstream ss;
+    ss << e;
+    Engine f;
+    ss >> f;
+    BOOST_TEST(e == f);
+}
+
+// Checks that seeding from a range with fewer values than the state throws.
+template <typename Engine>
+void check_short_seed()
+{
+    const std::size_t k = Engine().state().size();
+    const std::vector<typename Engine::result_type> v(k - 1, 1);
+    for (std::size_t n = 0; n < k; ++n)
+    {
+        Engine e;
+        BOOST_TEST_THROWS(e.seed(v.begin(), v.begin() + n), std::invalid_argument);
+        auto first = v.begin();
+        BOOST_TEST_THROWS(Engine(first, v.begin() + n), std::invalid_argument);
+    }
+}
+
+#endif // BOOST_RANDOM_TEST_COMP_MWCX_HPP
diff --git a/test/test_mwc128x.cpp b/test/test_mwc128x.cpp
new file mode 100644
index 000000000..a35778953
--- /dev/null
+++ b/test/test_mwc128x.cpp
@@ -0,0 +1,25 @@
+/* test_mwc128x.cpp
+ *
+ * Copyright Sebastiano Vigna 2026
+ * Distributed under the Boost Software License, Version 1.0. (See
+ * accompanying file LICENSE_1_0.txt or copy at
+ * http://www.boost.org/LICENSE_1_0.txt)
+ *
+ */
+
+#include <boost/random/mwcx.hpp>
+#include <cstdint>
+
+#define BOOST_RANDOM_URNG boost::random::mwc128x
+#define BOOST_RANDOM_CPP11_URNG
+#define BOOST_RANDOM_PROVIDED_SEED_TYPE std::uint64_t
+
+// Values computed with the reference implementation (https://github.com/vigna/MWCPRNG)
+// from the initial states produced by the default constructor and by seed_seq
+#define BOOST_RANDOM_VALIDATION_VALUE UINT64_C(0x672896149bed570a)
+#define BOOST_RANDOM_SEED_SEQ_VALIDATION_VALUE UINT64_C(0xab38168a1714dac3)
+
+// Larger than the distance below which discard() steps the generator
+#define BOOST_RANDOM_DISCARD_COUNT2 1000000
+
+#include "test_generator.ipp"
diff --git a/test/test_mwc128x32.cpp b/test/test_mwc128x32.cpp
new file mode 100644
index 000000000..4bd61a67f
--- /dev/null
+++ b/test/test_mwc128x32.cpp
@@ -0,0 +1,25 @@
+/* test_mwc128x32.cpp
+ *
+ * Copyright Sebastiano Vigna 2026
+ * Distributed under the Boost Software License, Version 1.0. (See
+ * accompanying file LICENSE_1_0.txt or copy at
+ * http://www.boost.org/LICENSE_1_0.txt)
+ *
+ */
+
+#include <boost/random/mwcx.hpp>
+#include <cstdint>
+
+#define BOOST_RANDOM_URNG boost::random::mwc128x32
+#define BOOST_RANDOM_CPP11_URNG
+#define BOOST_RANDOM_PROVIDED_SEED_TYPE std::uint64_t
+
+// Values computed with the reference implementation (https://github.com/vigna/MWCPRNG)
+// from the initial states produced by the default constructor and by seed_seq
+#define BOOST_RANDOM_VALIDATION_VALUE UINT32_C(0xce976888)
+#define BOOST_RANDOM_SEED_SEQ_VALIDATION_VALUE UINT32_C(0xcb40dfb2)
+
+// Larger than the distance below which discard() steps the generator
+#define BOOST_RANDOM_DISCARD_COUNT2 1000000
+
+#include "test_generator.ipp"
diff --git a/test/test_mwc192x.cpp b/test/test_mwc192x.cpp
new file mode 100644
index 000000000..4d2a46ba6
--- /dev/null
+++ b/test/test_mwc192x.cpp
@@ -0,0 +1,25 @@
+/* test_mwc192x.cpp
+ *
+ * Copyright Sebastiano Vigna 2026
+ * Distributed under the Boost Software License, Version 1.0. (See
+ * accompanying file LICENSE_1_0.txt or copy at
+ * http://www.boost.org/LICENSE_1_0.txt)
+ *
+ */
+
+#include <boost/random/mwcx.hpp>
+#include <cstdint>
+
+#define BOOST_RANDOM_URNG boost::random::mwc192x
+#define BOOST_RANDOM_CPP11_URNG
+#define BOOST_RANDOM_PROVIDED_SEED_TYPE std::uint64_t
+
+// Values computed with the reference implementation (https://github.com/vigna/MWCPRNG)
+// from the initial states produced by the default constructor and by seed_seq
+#define BOOST_RANDOM_VALIDATION_VALUE UINT64_C(0xf9849a66b8c3b458)
+#define BOOST_RANDOM_SEED_SEQ_VALIDATION_VALUE UINT64_C(0xb9674f3c536843dc)
+
+// Larger than the distance below which discard() steps the generator
+#define BOOST_RANDOM_DISCARD_COUNT2 1000000
+
+#include "test_generator.ipp"
diff --git a/test/test_mwc64x32.cpp b/test/test_mwc64x32.cpp
new file mode 100644
index 000000000..52d8ca240
--- /dev/null
+++ b/test/test_mwc64x32.cpp
@@ -0,0 +1,25 @@
+/* test_mwc64x32.cpp
+ *
+ * Copyright Sebastiano Vigna 2026
+ * Distributed under the Boost Software License, Version 1.0. (See
+ * accompanying file LICENSE_1_0.txt or copy at
+ * http://www.boost.org/LICENSE_1_0.txt)
+ *
+ */
+
+#include <boost/random/mwcx.hpp>
+#include <cstdint>
+
+#define BOOST_RANDOM_URNG boost::random::mwc64x32
+#define BOOST_RANDOM_CPP11_URNG
+#define BOOST_RANDOM_PROVIDED_SEED_TYPE std::uint64_t
+
+// Values computed with the reference implementation (https://github.com/vigna/MWCPRNG)
+// from the initial states produced by the default constructor and by seed_seq
+#define BOOST_RANDOM_VALIDATION_VALUE UINT32_C(0x2bd590a1)
+#define BOOST_RANDOM_SEED_SEQ_VALIDATION_VALUE UINT32_C(0xea06b0bd)
+
+// Larger than the distance below which discard() steps the generator
+#define BOOST_RANDOM_DISCARD_COUNT2 1000000
+
+#include "test_generator.ipp"
diff --git a/test/test_mwcx_mul_add.cpp b/test/test_mwcx_mul_add.cpp
new file mode 100644
index 000000000..20f09ad12
--- /dev/null
+++ b/test/test_mwcx_mul_add.cpp
@@ -0,0 +1,96 @@
+/*
+ * Copyright Sebastiano Vigna 2026.
+ * Distributed under the Boost Software License, Version 1.0. (See
+ * accompanying file LICENSE_1_0.txt or copy at
+ * http://www.boost.org/LICENSE_1_0.txt)
+ *
+ * Checks the 64 x 64 + 64 -> 128-bit operation of the scrambled
+ * multiply-with-carry generators against long multiplication on 16-bit
+ * limbs (and against unsigned __int128, if available).
+ */
+
+#include <boost/random/detail/mwcx_arithmetic.hpp>
+#include <boost/core/lightweight_test.hpp>
+#include <cstdint>
+
+using std::uint64_t;
+
+// a * x + c by long multiplication on 16-bit limbs
+static void reference(const uint64_t a, const uint64_t x, const uint64_t c, uint64_t& lo, uint64_t& hi)
+{
+    uint64_t r[8] = { 0 };
+    for (int i = 0; i < 4; ++i)
+    {
+        r[i] = (c >> (16 * i)) & 0xffff;
+    }
+
+    for (int i = 0; i < 4; ++i)
+    {
+        uint64_t carry = 0;
+        for (int j = 0; j < 4; ++j)
+        {
+            const uint64_t t = ((a >> (16 * i)) & 0xffff) * ((x >> (16 * j)) & 0xffff) + r[i + j] + carry;
+            r[i + j] = t & 0xffff;
+            carry = t >> 16;
+        }
+
+        for (int k = i + 4; carry != 0; ++k)
+        {
+            const uint64_t t = r[k] + carry;
+            r[k] = t & 0xffff;
+            carry = t >> 16;
+        }
+    }
+
+    lo = r[0] | r[1] << 16 | r[2] << 32 | r[3] << 48;
+    hi = r[4] | r[5] << 16 | r[6] << 32 | r[7] << 48;
+}
+
+static void check(const uint64_t a, const uint64_t x, const uint64_t c)
+{
+    uint64_t hi;
+    const uint64_t lo = boost::random::detail::mwcx_mul_add(a, x, c, hi);
+    uint64_t rlo, rhi;
+    reference(a, x, c, rlo, rhi);
+    BOOST_TEST_EQ(lo, rlo);
+    BOOST_TEST_EQ(hi, rhi);
+#if defined(__SIZEOF_INT128__)
+    __extension__ typedef unsigned __int128 uint128_type;
+    const uint128_type t = static_cast<uint128_type>(a) * x + c;
+    BOOST_TEST_EQ(lo, static_cast<uint64_t>(t));
+    BOOST_TEST_EQ(hi, static_cast<uint64_t>(t >> 64));
+#endif
+}
+
+int main()
+{
+    const uint64_t max = UINT64_C(0xffffffffffffffff);
+    const uint64_t edge[] = { 0, 1, 2, UINT64_C(0xffffffff), UINT64_C(0x100000000), UINT64_C(0x8000000000000000), max - 1, max };
+    for (const uint64_t a : edge)
+    {
+        for (const uint64_t x : edge)
+        {
+            for (const uint64_t c : edge)
+            {
+                check(a, x, c);
+            }
+        }
+    }
+
+    uint64_t s = UINT64_C(0x9e3779b97f4a7c15);
+    for (int i = 0; i < 1000000; ++i)
+    {
+        s ^= s << 13;
+        s ^= s >> 7;
+        s ^= s << 17;
+        const uint64_t a = s * UINT64_C(0x2545F4914F6CDD1D);
+        const uint64_t x = s;
+        const uint64_t c = s * UINT64_C(0x9E3779B97F4A7C15);
+        check(a, x, c);
+        // The multipliers of the generators
+        check(UINT64_C(0xff4e24dd418f5433), x, c);
+        check(UINT64_C(0xff6da236274e252e), x, c);
+    }
+
+    return boost::report_errors();
+}
diff --git a/test/test_mwcx_mul_add_portable.cpp b/test/test_mwcx_mul_add_portable.cpp
new file mode 100644
index 000000000..722d708c7
--- /dev/null
+++ b/test/test_mwcx_mul_add_portable.cpp
@@ -0,0 +1,11 @@
+/*
+ * Copyright Sebastiano Vigna 2026.
+ * Distributed under the Boost Software License, Version 1.0. (See
+ * accompanying file LICENSE_1_0.txt or copy at
+ * http://www.boost.org/LICENSE_1_0.txt)
+ *
+ * The tests of test_mwcx_mul_add.cpp with the portable 64 x 64 -> 128-bit multiplication.
+ */
+
+#define BOOST_RANDOM_MWCX_PORTABLE_MUL
+#include "test_mwcx_mul_add.cpp"