diff --git a/math/mathcore/inc/Math/RanluxppEngine.h b/math/mathcore/inc/Math/RanluxppEngine.h
index bfef96bc7735c..d790f75027988 100644
--- a/math/mathcore/inc/Math/RanluxppEngine.h
+++ b/math/mathcore/inc/Math/RanluxppEngine.h
@@ -2,7 +2,7 @@
 // Author: Jonas Hahnfeld 11/2020
 
 /*************************************************************************
- * Copyright (C) 1995-2020, Rene Brun and Fons Rademakers.               *
+ * Copyright (C) 1995-2021, Rene Brun and Fons Rademakers.               *
  * All rights reserved.                                                  *
  *                                                                       *
  * For the licensing terms see $ROOTSYS/LICENSE.                         *
@@ -20,14 +20,15 @@
 namespace ROOT {
 namespace Math {
 
-template <int w, int p>
+template <int w, int p, int u = 0>
 class RanluxppEngineImpl;
 
 template <int p>
 class RanluxppEngine final : public TRandomEngine {
 
 private:
-   std::unique_ptr<RanluxppEngineImpl<48, p>> fImpl;
+   using ImplType = RanluxppEngineImpl<48, p>;
+   std::unique_ptr<ImplType> fImpl;
 
 public:
    RanluxppEngine(uint64_t seed = 314159265);
@@ -55,6 +56,246 @@ using RanluxppEngine2048 = RanluxppEngine<2048>;
 extern template class RanluxppEngine<24>;
 extern template class RanluxppEngine<2048>;
 
+
+template <int p>
+class RanluxppCompatEngineJames final : public TRandomEngine {
+
+private:
+   using ImplType = RanluxppEngineImpl<24, p>;
+   std::unique_ptr<ImplType> fImpl;
+
+public:
+   RanluxppCompatEngineJames(uint64_t seed = 314159265);
+   virtual ~RanluxppCompatEngineJames();
+
+   /// Generate a floating point random number with 24 bits of randomness
+   double Rndm() override;
+   /// Generate a floating point random number (non-virtual method)
+   double operator()();
+   /// Generate a random integer value with 24 bits
+   uint64_t IntRndm();
+
+   /// Initialize and seed the state of the generator
+   void SetSeed(uint64_t seed);
+   /// Skip `n` random numbers without generating them
+   void Skip(uint64_t n);
+
+   /// Get name of the generator
+   static const char *Name() { return "RanluxppCompatJames"; }
+};
+
+/// Compatibility engine for original RANLUX implementation by James, luxury
+/// level 3 (p = 223). The sequence of numbers also matches `gsl_rng_ranlux`.
+using RanluxppCompatEngineJamesP3 = RanluxppCompatEngineJames<223>;
+/// Compatibility engine for original RANLUX implementation by James, luxury
+/// level 4 (p = 389). The sequence of numbers also matches `gsl_rng_ranlux389`.
+using RanluxppCompatEngineJamesP4 = RanluxppCompatEngineJames<389>;
+
+extern template class RanluxppCompatEngineJames<223>;
+extern template class RanluxppCompatEngineJames<389>;
+
+
+/// Compatibility engine for `gsl_rng_ranlxs*` from the GNU Scientific Library.
+template <int p>
+class RanluxppCompatEngineGslRanlxs final : public TRandomEngine {
+
+private:
+   using ImplType = RanluxppEngineImpl<24, p>;
+   std::unique_ptr<ImplType> fImpl;
+
+public:
+   RanluxppCompatEngineGslRanlxs(uint64_t seed = 1);
+   virtual ~RanluxppCompatEngineGslRanlxs();
+
+   /// Generate a floating point random number with 24 bits of randomness
+   double Rndm() override;
+   /// Generate a floating point random number (non-virtual method)
+   double operator()();
+   /// Generate a random integer value with 24 bits
+   uint64_t IntRndm();
+
+   /// Initialize and seed the state of the generator
+   void SetSeed(uint64_t seed);
+   /// Skip `n` random numbers without generating them
+   void Skip(uint64_t n);
+
+   /// Get name of the generator
+   static const char *Name() { return "RanluxppCompatGslRanlxs"; }
+};
+
+using RanluxppCompatEngineGslRanlxs0 = RanluxppCompatEngineGslRanlxs<218>;
+using RanluxppCompatEngineGslRanlxs1 = RanluxppCompatEngineGslRanlxs<404>;
+using RanluxppCompatEngineGslRanlxs2 = RanluxppCompatEngineGslRanlxs<794>;
+
+extern template class RanluxppCompatEngineGslRanlxs<218>;
+extern template class RanluxppCompatEngineGslRanlxs<404>;
+extern template class RanluxppCompatEngineGslRanlxs<794>;
+
+
+/// Compatibility engine for `gsl_rng_ranlxd*` from the GNU Scientific Library.
+template <int p>
+class RanluxppCompatEngineGslRanlxd final : public TRandomEngine {
+
+private:
+   using ImplType = RanluxppEngineImpl<48, p>;
+   std::unique_ptr<ImplType> fImpl;
+
+public:
+   RanluxppCompatEngineGslRanlxd(uint64_t seed = 1);
+   virtual ~RanluxppCompatEngineGslRanlxd();
+
+   /// Generate a floating point random number with 48 bits of randomness
+   double Rndm() override;
+   /// Generate a floating point random number (non-virtual method)
+   double operator()();
+   /// Generate a random integer value with 48 bits
+   uint64_t IntRndm();
+
+   /// Initialize and seed the state of the generator
+   void SetSeed(uint64_t seed);
+   /// Skip `n` random numbers without generating them
+   void Skip(uint64_t n);
+
+   /// Get name of the generator
+   static const char *Name() { return "RanluxppCompatGslRanlxd"; }
+};
+
+using RanluxppCompatEngineGslRanlxd1 = RanluxppCompatEngineGslRanlxd<404>;
+using RanluxppCompatEngineGslRanlxd2 = RanluxppCompatEngineGslRanlxd<794>;
+
+extern template class RanluxppCompatEngineGslRanlxd<404>;
+extern template class RanluxppCompatEngineGslRanlxd<794>;
+
+
+template <int w, int p>
+class RanluxppCompatEngineLuescherImpl;
+
+/// Compatibility engine for Lüscher's ranlxs implementation written in C.
+template <int p>
+class RanluxppCompatEngineLuescherRanlxs final : public TRandomEngine {
+
+private:
+   using ImplType = RanluxppCompatEngineLuescherImpl<24, p>;
+   std::unique_ptr<ImplType> fImpl;
+
+public:
+   RanluxppCompatEngineLuescherRanlxs(uint64_t seed = 314159265);
+   virtual ~RanluxppCompatEngineLuescherRanlxs();
+
+   /// Generate a floating point random number with 24 bits of randomness
+   double Rndm() override;
+   /// Generate a floating point random number (non-virtual method)
+   double operator()();
+   /// Generate a random integer value with 24 bits
+   uint64_t IntRndm();
+
+   /// Initialize and seed the state of the generator
+   void SetSeed(uint64_t seed);
+   /// Skip `n` random numbers without generating them
+   void Skip(uint64_t n);
+
+   /// Get name of the generator
+   static const char *Name() { return "RanluxppCompatLuescherRanlxs"; }
+};
+
+using RanluxppCompatEngineLuescherRanlxs0 = RanluxppCompatEngineLuescherRanlxs<218>;
+using RanluxppCompatEngineLuescherRanlxs1 = RanluxppCompatEngineLuescherRanlxs<404>;
+using RanluxppCompatEngineLuescherRanlxs2 = RanluxppCompatEngineLuescherRanlxs<794>;
+
+extern template class RanluxppCompatEngineLuescherRanlxs<218>;
+extern template class RanluxppCompatEngineLuescherRanlxs<404>;
+extern template class RanluxppCompatEngineLuescherRanlxs<794>;
+
+
+/// Compatibility engine for Lüscher's ranlxd implementation written in C.
+template <int p>
+class RanluxppCompatEngineLuescherRanlxd final : public TRandomEngine {
+
+private:
+   using ImplType = RanluxppCompatEngineLuescherImpl<48, p>;
+   std::unique_ptr<ImplType> fImpl;
+
+public:
+   RanluxppCompatEngineLuescherRanlxd(uint64_t seed = 314159265);
+   virtual ~RanluxppCompatEngineLuescherRanlxd();
+
+   /// Generate a floating point random number with 48 bits of randomness
+   double Rndm() override;
+   /// Generate a floating point random number (non-virtual method)
+   double operator()();
+   /// Generate a random integer value with 48 bits
+   uint64_t IntRndm();
+
+   /// Initialize and seed the state of the generator
+   void SetSeed(uint64_t seed);
+   /// Skip `n` random numbers without generating them
+   void Skip(uint64_t n);
+
+   /// Get name of the generator
+   static const char *Name() { return "RanluxppCompatLuescherRanlxd"; }
+};
+using RanluxppCompatEngineLuescherRanlxd1 = RanluxppCompatEngineLuescherRanlxd<404>;
+using RanluxppCompatEngineLuescherRanlxd2 = RanluxppCompatEngineLuescherRanlxd<794>;
+
+extern template class RanluxppCompatEngineLuescherRanlxd<404>;
+extern template class RanluxppCompatEngineLuescherRanlxd<794>;
+
+
+/// Compatibility engine for `std::ranlux24` from the C++ standard.
+class RanluxppCompatEngineStdRanlux24 final : public TRandomEngine {
+
+private:
+   using ImplType = RanluxppEngineImpl<24, 223, 23>;
+   std::unique_ptr<ImplType> fImpl;
+
+public:
+   RanluxppCompatEngineStdRanlux24(uint64_t seed = 19780503);
+   virtual ~RanluxppCompatEngineStdRanlux24();
+
+   /// Generate a floating point random number with 24 bits of randomness
+   double Rndm() override;
+   /// Generate a floating point random number (non-virtual method)
+   double operator()();
+   /// Generate a random integer value with 24 bits
+   uint64_t IntRndm();
+
+   /// Initialize and seed the state of the generator
+   void SetSeed(uint64_t seed);
+   /// Skip `n` random numbers without generating them
+   void Skip(uint64_t n);
+
+   /// Get name of the generator
+   static const char *Name() { return "RanluxppCompatStdRanlux24"; }
+};
+
+
+/// Compatibility engine for `std::ranlux48` from the C++ standard.
+class RanluxppCompatEngineStdRanlux48 final : public TRandomEngine {
+
+private:
+   using ImplType = RanluxppEngineImpl<48, 2 * 389, 11>;
+   std::unique_ptr<ImplType> fImpl;
+
+public:
+   RanluxppCompatEngineStdRanlux48(uint64_t seed = 19780503);
+   virtual ~RanluxppCompatEngineStdRanlux48();
+
+   /// Generate a floating point random number with 48 bits of randomness
+   double Rndm() override;
+   /// Generate a floating point random number (non-virtual method)
+   double operator()();
+   /// Generate a random integer value with 48 bits
+   uint64_t IntRndm();
+
+   /// Initialize and seed the state of the generator
+   void SetSeed(uint64_t seed);
+   /// Skip `n` random numbers without generating them
+   void Skip(uint64_t n);
+
+   /// Get name of the generator
+   static const char *Name() { return "RanluxppCompatStdRanlux48"; }
+};
+
 } // end namespace Math
 } // end namespace ROOT
 
diff --git a/math/mathcore/src/RanluxppEngineImpl.cxx b/math/mathcore/src/RanluxppEngineImpl.cxx
index 83ecbd7f576d3..100f8d8638765 100644
--- a/math/mathcore/src/RanluxppEngineImpl.cxx
+++ b/math/mathcore/src/RanluxppEngineImpl.cxx
@@ -2,7 +2,7 @@
 // Author: Jonas Hahnfeld 11/2020
 
 /*************************************************************************
- * Copyright (C) 1995-2020, Rene Brun and Fons Rademakers.               *
+ * Copyright (C) 1995-2021, Rene Brun and Fons Rademakers.               *
  * All rights reserved.                                                  *
  *                                                                       *
  * For the licensing terms see $ROOTSYS/LICENSE.                         *
@@ -49,6 +49,18 @@ namespace {
 // Variable templates are a feature of C++14, use the older technique of having
 // a static member in a template class.
 
+// The coefficients have been determined using Python, and in parts compared to the values given by Sibidanov.
+//
+//     >>> def print_hex(a):
+//     ...     while a > 0:
+//     ...         print('{0:#018x}'.format(a & 0xffffffffffffffff))
+//     ...         a >>= 64
+//     ...
+//     >>> m = 2 ** 576 - 2 ** 240 + 1
+//     >>> a = m - (m - 1) // 2 ** 24
+//     >>> kA = pow(a, <w>, m)
+//     >>> print_hex(kA)
+
 template <int p>
 struct RanluxppData;
 
@@ -56,15 +68,73 @@ template <>
 struct RanluxppData<24> {
    static const uint64_t kA[9];
 };
+// Also given by Sibidanov
 const uint64_t RanluxppData<24>::kA[] = {
    0x0000000000000000, 0x0000000000000000, 0x0000000000010000, 0xfffe000000000000, 0xffffffffffffffff,
    0xffffffffffffffff, 0xffffffffffffffff, 0xfffffffeffffffff, 0xffffffffffffffff,
 };
 
+template <>
+struct RanluxppData<218> {
+   static const uint64_t kA[9];
+};
+const uint64_t RanluxppData<218>::kA[] = {
+	0xf445fffffffffd94, 0xfffffd74ffffffff, 0x000000000ba5ffff, 0xfc76000000000942, 0xfffffaaaffffffff,
+	0x0000000000b0ffff, 0x027b0000000007d1, 0xfffff96000000000, 0xfffffffff8e4ffff,
+};
+
+template <>
+struct RanluxppData<223> {
+   static const uint64_t kA[9];
+};
+// Also given by Sibidanov
+const uint64_t RanluxppData<223>::kA[] = {
+   0x0000000ba6000000, 0x0a00000000094200, 0xffeef0fffffffffa, 0xfffffffe25ffffff, 0x7b0000000007d0ff,
+   0xfff9600000000002, 0xfffffff8e4ffffff, 0xba00000000026cff, 0x00028b000000000b,
+};
+
+template <>
+struct RanluxppData<389> {
+   static const uint64_t kA[9];
+};
+// Also given by Sibidanov
+const uint64_t RanluxppData<389>::kA[] = {
+   0x00002ecac9000000, 0x740000002c389600, 0xb9c8a6ffffffe525, 0xfffff593cfffffff, 0xab0000001e93f2ff,
+   0xe4ab160000000d92, 0xffffdf6604ffffff, 0x020000000b9242ff, 0x0df0600000002ee0,
+};
+
+template <>
+struct RanluxppData<404> {
+   static const uint64_t kA[9];
+};
+const uint64_t RanluxppData<404>::kA[] = {
+	0x2eabffffffc9d08b, 0x00012612ffffff99, 0x0000007c3ebe0000, 0x353600000047bba1, 0xffd3c769ffffffd1,
+	0x0000001ada8bffff, 0x6c30000000463759, 0xffb2a1440000000a, 0xffffffc634beffff,
+};
+
+template <>
+struct RanluxppData<778> {
+   static const uint64_t kA[9];
+};
+const uint64_t RanluxppData<778>::kA[] = {
+   0x872de42d9dca512b, 0xdbf015ea1662f8a0, 0x01f48f0d28482e96, 0x392fca0b3be2ae04, 0xed00881af896ce54,
+   0x14f0a768664013f3, 0x9489f52deb1f7f80, 0x72139804e09c0f37, 0x2146b0bb92a2f9a4,
+};
+
+template <>
+struct RanluxppData<794> {
+   static const uint64_t kA[9];
+};
+const uint64_t RanluxppData<794>::kA[] = {
+	0x428df7227a2ca7c9, 0xde32225faaa74b1a, 0x4b9d965ca1ebd668, 0x78d15f59e58e2aff, 0x240fea15e99d075f,
+	0xfe0b70f2d7b7d169, 0x75a535f4c41d51fb, 0x1a5ef0b7233b93e1, 0xbc787ca783d5d5a9,
+};
+
 template <>
 struct RanluxppData<2048> {
    static const uint64_t kA[9];
 };
+// Also given by Sibidanov
 const uint64_t RanluxppData<2048>::kA[] = {
    0xed7faa90747aaad9, 0x4cec2c78af55c101, 0xe64dcb31c48228ec, 0x6d8a15a13bee7cb0, 0x20b2ca60cb78c509,
    0x256c3d3c662ea36c, 0xff74e54107684ed2, 0x492edfcc0cc8e753, 0xb48c187cf5b22097,
@@ -75,8 +145,10 @@ const uint64_t RanluxppData<2048>::kA[] = {
 namespace ROOT {
 namespace Math {
 
-template <int w, int p>
+template <int w, int p, int u>
 class RanluxppEngineImpl {
+   // Needs direct access to private members to initialize its four states.
+   friend class RanluxppCompatEngineLuescherImpl<w, p>;
 
 private:
    uint64_t fState[9]; ///< RANLUX state of the generator
@@ -84,18 +156,31 @@ class RanluxppEngineImpl {
    int fPosition = 0;  ///< Current position in bits
 
    static constexpr const uint64_t *kA = RanluxppData<p>::kA;
-   static constexpr int kMaxPos = 9 * 64;
+   static constexpr int kMaxPos = (u == 0) ? 9 * 64 : u * w;
+   static_assert(kMaxPos <= 576, "maximum position larger than 576 bits");
 
-   /// Produce next block of random bits
-   void Advance()
+   /// Advance with given multiplier
+   void Advance(const uint64_t *a)
    {
       uint64_t lcg[9];
       to_lcg(fState, fCarry, lcg);
-      mulmod(kA, lcg);
+      mulmod(a, lcg);
       to_ranlux(lcg, fState, fCarry);
       fPosition = 0;
    }
 
+   /// Produce next block of random bits
+   void Advance()
+   {
+      Advance(kA);
+   }
+
+   /// Skip 24 RANLUX numbers
+   void Skip24()
+   {
+      Advance(RanluxppData<24>::kA);
+   }
+
 public:
    /// Return the next random bits, generate a new block if necessary
    uint64_t NextRandomBits()
@@ -128,8 +213,84 @@ class RanluxppEngineImpl {
       return bits * div;
    }
 
-   /// Initialize and seed the state of the generator
-   void SetSeed(uint64_t s)
+   /// Initialize and seed the state of the generator as in James' implementation
+   void SetSeedJames(uint64_t s)
+   {
+      // Multiplicative Congruential generator using formula constants of L'Ecuyer
+      // as described in "A review of pseudorandom number generators" (Fred James)
+      // published in Computer Physics Communications 60 (1990) pages 329-344.
+      int64_t seed = s;
+      auto next = [&]() {
+         const int a = 0xd1a4, b = 0x9c4e, c = 0x2fb3, d = 0x7fffffab;
+         int64_t k = seed / a;
+         seed = b * (seed - k * a) - k * c ;
+         if (seed < 0) seed += d;
+         return seed & 0xffffff;
+      };
+
+      // Iteration is reversed because the first number from the MCG goes to the
+      // highest position.
+      for (int i = 6; i >= 0; i -= 3) {
+         uint64_t r[8];
+         for (int j = 0; j < 8; j++) {
+            r[j] = next();
+         }
+
+         fState[i+0] = r[7] + (r[6] << 24) + (r[5] << 48);
+         fState[i+1] = (r[5] >> 16) + (r[4] << 8) + (r[3] << 32) + (r[2] << 56);
+         fState[i+2] = (r[2] >> 8) + (r[1] << 16) + (r[0] << 40);
+      }
+      fCarry = !seed;
+
+      Skip24();
+   }
+
+   /// Initialize and seed the state of the generator as in gsl_rng_ranlx*
+   void SetSeedGsl(uint32_t s, bool ranlxd)
+   {
+      if (s == 0) {
+         // The default seed for gsl_rng_ranlx* is 1.
+         s = 1;
+      }
+
+      uint32_t bits = s;
+      auto next_bit = [&]() {
+         int b13 = (bits >> 18) & 0x1;
+         int b31 = bits & 0x1;
+         uint32_t bn = b13 ^ b31;
+         bits = (bn << 30) + (bits >> 1);
+         return b31;
+      };
+      auto next = [&]() {
+         uint64_t ix = 0;
+         for (int i = 0; i < 48; i++) {
+            int iy = next_bit();
+            if (ranlxd) {
+               iy = (iy + 1) % 2;
+            }
+            ix = 2 * ix + iy;
+         }
+         return ix;
+      };
+
+      for (int i = 0; i < 9; i += 3) {
+         uint64_t r[4];
+         for (int j = 0; j < 4; j++) {
+            r[j] = next();
+         }
+
+         fState[i+0] = r[0] + (r[1] << 48);
+         fState[i+1] = (r[1] >> 16) + (r[2] << 32);
+         fState[i+2] = (r[2] >> 32) + (r[3] << 16);
+      }
+
+      fCarry = 0;
+      fPosition = 0;
+      Advance();
+   }
+
+   /// Initialize and seed the state of the generator as proposed by Sibidanov
+   void SetSeedSibidanov(uint64_t s)
    {
       uint64_t lcg[9];
       lcg[0] = 1;
@@ -149,6 +310,61 @@ class RanluxppEngineImpl {
       fPosition = 0;
    }
 
+   /// Initialize and seed the state of the generator as described by the C++ standard
+   void SetSeedStd24(uint64_t s)
+   {
+      // Seed LCG with given parameters.
+      uint64_t seed = s;
+      const uint64_t a = 40014, m = 2147483563;
+      auto next = [&]() {
+         seed = (a * seed) % m;
+         return seed & 0xffffff;
+      };
+
+      for (int i = 0; i < 9; i += 3) {
+         uint64_t r[8];
+         for (int j = 0; j < 8; j++) {
+            r[j] = next();
+         }
+
+         fState[i+0] = r[0] + (r[1] << 24) + (r[2] << 48);
+         fState[i+1] = (r[2] >> 16) + (r[3] << 8) + (r[4] << 32) + (r[5] << 56);
+         fState[i+2] = (r[5] >> 8) + (r[6] << 16) + (r[7] << 40);
+      }
+      fCarry = !seed;
+
+      Skip24();
+   }
+
+   /// Initialize and seed the state of the generator as described by the C++ standard
+   void SetSeedStd48(uint64_t s)
+   {
+      // Seed LCG with given parameters.
+      uint64_t seed = s;
+      const uint64_t a = 40014, m = 2147483563;
+      auto next = [&]() {
+         seed = (a * seed) % m;
+         uint64_t result = seed;
+         seed = (a * seed) % m;
+         result += seed << 32;
+         return result & 0xffffffffffff;
+      };
+
+      for (int i = 0; i < 9; i += 3) {
+         uint64_t r[4];
+         for (int j = 0; j < 4; j++) {
+            r[j] = next();
+         }
+
+         fState[i+0] = r[0] + (r[1] << 48);
+         fState[i+1] = (r[1] >> 16) + (r[2] << 32);
+         fState[i+2] = (r[2] >> 32) + (r[3] << 16);
+      }
+      fCarry = !seed;
+
+      Skip24();
+   }
+
    /// Skip `n` random numbers without generating them
    void Skip(uint64_t n)
    {
@@ -183,9 +399,9 @@ class RanluxppEngineImpl {
 };
 
 template <int p>
-RanluxppEngine<p>::RanluxppEngine(uint64_t seed) : fImpl(new RanluxppEngineImpl<48, p>)
+RanluxppEngine<p>::RanluxppEngine(uint64_t seed) : fImpl(new ImplType)
 {
-   fImpl->SetSeed(seed);
+   this->SetSeed(seed);
 }
 
 template <int p>
@@ -212,7 +428,7 @@ uint64_t RanluxppEngine<p>::IntRndm()
 template <int p>
 void RanluxppEngine<p>::SetSeed(uint64_t seed)
 {
-   fImpl->SetSeed(seed);
+   fImpl->SetSeedSibidanov(seed);
 }
 
 template <int p>
@@ -224,5 +440,376 @@ void RanluxppEngine<p>::Skip(uint64_t n)
 template class RanluxppEngine<24>;
 template class RanluxppEngine<2048>;
 
+
+template <int p>
+RanluxppCompatEngineJames<p>::RanluxppCompatEngineJames(uint64_t seed) : fImpl(new ImplType)
+{
+   this->SetSeed(seed);
+}
+
+template <int p>
+RanluxppCompatEngineJames<p>::~RanluxppCompatEngineJames() = default;
+
+template <int p>
+double RanluxppCompatEngineJames<p>::Rndm()
+{
+   return (*this)();
+}
+
+template <int p>
+double RanluxppCompatEngineJames<p>::operator()()
+{
+   return fImpl->NextRandomFloat();
+}
+
+template <int p>
+uint64_t RanluxppCompatEngineJames<p>::IntRndm()
+{
+   return fImpl->NextRandomBits();
+}
+
+template <int p>
+void RanluxppCompatEngineJames<p>::SetSeed(uint64_t seed)
+{
+   fImpl->SetSeedJames(seed);
+}
+
+template <int p>
+void RanluxppCompatEngineJames<p>::Skip(uint64_t n)
+{
+   fImpl->Skip(n);
+}
+
+template class RanluxppCompatEngineJames<223>;
+template class RanluxppCompatEngineJames<389>;
+
+
+template <int p>
+RanluxppCompatEngineGslRanlxs<p>::RanluxppCompatEngineGslRanlxs(uint64_t seed) : fImpl(new ImplType)
+{
+   this->SetSeed(seed);
+}
+
+template <int p>
+RanluxppCompatEngineGslRanlxs<p>::~RanluxppCompatEngineGslRanlxs() = default;
+
+template <int p>
+double RanluxppCompatEngineGslRanlxs<p>::Rndm()
+{
+   return (*this)();
+}
+
+template <int p>
+double RanluxppCompatEngineGslRanlxs<p>::operator()()
+{
+   return fImpl->NextRandomFloat();
+}
+
+template <int p>
+uint64_t RanluxppCompatEngineGslRanlxs<p>::IntRndm()
+{
+   return fImpl->NextRandomBits();
+}
+
+template <int p>
+void RanluxppCompatEngineGslRanlxs<p>::SetSeed(uint64_t seed)
+{
+   fImpl->SetSeedGsl(seed, /*ranlxd=*/false);
+}
+
+template <int p>
+void RanluxppCompatEngineGslRanlxs<p>::Skip(uint64_t n)
+{
+   fImpl->Skip(n);
+}
+
+template class RanluxppCompatEngineGslRanlxs<218>;
+template class RanluxppCompatEngineGslRanlxs<404>;
+template class RanluxppCompatEngineGslRanlxs<794>;
+
+
+template <int p>
+RanluxppCompatEngineGslRanlxd<p>::RanluxppCompatEngineGslRanlxd(uint64_t seed) : fImpl(new ImplType)
+{
+   this->SetSeed(seed);
+}
+
+template <int p>
+RanluxppCompatEngineGslRanlxd<p>::~RanluxppCompatEngineGslRanlxd() = default;
+
+template <int p>
+double RanluxppCompatEngineGslRanlxd<p>::Rndm()
+{
+   return (*this)();
+}
+
+template <int p>
+double RanluxppCompatEngineGslRanlxd<p>::operator()()
+{
+   return fImpl->NextRandomFloat();
+}
+
+template <int p>
+uint64_t RanluxppCompatEngineGslRanlxd<p>::IntRndm()
+{
+   return fImpl->NextRandomBits();
+}
+
+template <int p>
+void RanluxppCompatEngineGslRanlxd<p>::SetSeed(uint64_t seed)
+{
+   fImpl->SetSeedGsl(seed, /*ranlxd=*/true);
+}
+
+template <int p>
+void RanluxppCompatEngineGslRanlxd<p>::Skip(uint64_t n)
+{
+   fImpl->Skip(n);
+}
+
+template class RanluxppCompatEngineGslRanlxd<404>;
+template class RanluxppCompatEngineGslRanlxd<794>;
+
+
+template <int w, int p>
+class RanluxppCompatEngineLuescherImpl {
+
+private:
+  RanluxppEngineImpl<w, p> fStates[4]; ///< The states of this generator
+  int fNextState = 0;                  ///< The index of the next state
+
+public:
+   /// Return the next random bits, generate a new block if necessary
+   uint64_t NextRandomBits()
+   {
+      uint64_t bits = fStates[fNextState].NextRandomBits();
+      fNextState = (fNextState + 1) % 4;
+      return bits;
+   }
+
+   /// Return a floating point number, converted from the next random bits.
+   double NextRandomFloat()
+   {
+      double number = fStates[fNextState].NextRandomFloat();
+      fNextState = (fNextState + 1) % 4;
+      return number;
+   }
+
+   /// Initialize and seed the state of the generator as in Lüscher's ranlxs
+   void SetSeed(uint32_t s, bool ranlxd)
+   {
+      uint32_t bits = s;
+      auto next_bit = [&]() {
+         int b13 = (bits >> 18) & 0x1;
+         int b31 = bits & 0x1;
+         uint32_t bn = b13 ^ b31;
+         bits = (bn << 30) + (bits >> 1);
+         return b31;
+      };
+      auto next = [&]() {
+         uint64_t ix = 0;
+         for (int l = 0; l < 24; l++) {
+            ix = 2 * ix + next_bit();
+         }
+         return ix;
+      };
+
+      for (int i = 0; i < 4; i++) {
+         auto &state = fStates[i];
+         for (int j = 0; j < 9; j += 3) {
+            uint64_t r[8];
+            for (int m = 0; m < 8; m++) {
+               uint64_t ix = next();
+               // Lüscher's implementation uses k = (j / 3) * 8 + m, so only
+               // the value of m is important for (k % 4).
+               if ((!ranlxd && (m % 4) == i) || (ranlxd && (m % 4) != i)) {
+                  ix = 16777215 - ix;
+               }
+               r[m] = ix;
+            }
+
+            state.fState[j+0] = r[0] + (r[1] << 24) + (r[2] << 48);
+            state.fState[j+1] = (r[2] >> 16) + (r[3] << 8) + (r[4] << 32) + (r[5] << 56);
+            state.fState[j+2] = (r[5] >> 8) + (r[6] << 16) + (r[7] << 40);
+         }
+
+         state.fCarry = 0;
+         state.fPosition = 0;
+         state.Advance();
+      }
+
+      fNextState = 0;
+   }
+
+   /// Skip `n` random numbers without generating them
+   void Skip(uint64_t n)
+   {
+      uint64_t nPerState = n / 4;
+      int remainder = n % 4;
+      for (int i = 0; i < 4; i++) {
+         int idx = (fNextState + i) % 4;
+         uint64_t nForThisState = nPerState;
+         if (i < remainder) {
+            nForThisState++;
+         }
+         fStates[idx].Skip(nForThisState);
+      }
+      // Switch the next state according to the remainder.
+      fNextState = (fNextState + remainder) % 4;
+   }
+};
+
+template <int p>
+RanluxppCompatEngineLuescherRanlxs<p>::RanluxppCompatEngineLuescherRanlxs(uint64_t seed) : fImpl(new ImplType)
+{
+   this->SetSeed(seed);
+}
+
+template <int p>
+RanluxppCompatEngineLuescherRanlxs<p>::~RanluxppCompatEngineLuescherRanlxs() = default;
+
+template <int p>
+double RanluxppCompatEngineLuescherRanlxs<p>::Rndm()
+{
+   return (*this)();
+}
+
+template <int p>
+double RanluxppCompatEngineLuescherRanlxs<p>::operator()()
+{
+   return fImpl->NextRandomFloat();
+}
+
+template <int p>
+uint64_t RanluxppCompatEngineLuescherRanlxs<p>::IntRndm()
+{
+   return fImpl->NextRandomBits();
+}
+
+template <int p>
+void RanluxppCompatEngineLuescherRanlxs<p>::SetSeed(uint64_t seed)
+{
+   fImpl->SetSeed(seed, /*ranlxd=*/false);
+}
+
+template <int p>
+void RanluxppCompatEngineLuescherRanlxs<p>::Skip(uint64_t n)
+{
+   fImpl->Skip(n);
+}
+
+template class RanluxppCompatEngineLuescherRanlxs<218>;
+template class RanluxppCompatEngineLuescherRanlxs<404>;
+template class RanluxppCompatEngineLuescherRanlxs<794>;
+
+
+template <int p>
+RanluxppCompatEngineLuescherRanlxd<p>::RanluxppCompatEngineLuescherRanlxd(uint64_t seed) : fImpl(new ImplType)
+{
+   this->SetSeed(seed);
+}
+
+template <int p>
+RanluxppCompatEngineLuescherRanlxd<p>::~RanluxppCompatEngineLuescherRanlxd() = default;
+
+template <int p>
+double RanluxppCompatEngineLuescherRanlxd<p>::Rndm()
+{
+   return (*this)();
+}
+
+template <int p>
+double RanluxppCompatEngineLuescherRanlxd<p>::operator()()
+{
+   return fImpl->NextRandomFloat();
+}
+
+template <int p>
+uint64_t RanluxppCompatEngineLuescherRanlxd<p>::IntRndm()
+{
+   return fImpl->NextRandomBits();
+}
+
+template <int p>
+void RanluxppCompatEngineLuescherRanlxd<p>::SetSeed(uint64_t seed)
+{
+   fImpl->SetSeed(seed, /*ranlxd=*/true);
+}
+
+template <int p>
+void RanluxppCompatEngineLuescherRanlxd<p>::Skip(uint64_t n)
+{
+   fImpl->Skip(n);
+}
+
+template class RanluxppCompatEngineLuescherRanlxd<404>;
+template class RanluxppCompatEngineLuescherRanlxd<794>;
+
+
+RanluxppCompatEngineStdRanlux24::RanluxppCompatEngineStdRanlux24(uint64_t seed) : fImpl(new ImplType)
+{
+   this->SetSeed(seed);
+}
+
+RanluxppCompatEngineStdRanlux24::~RanluxppCompatEngineStdRanlux24() = default;
+
+double RanluxppCompatEngineStdRanlux24::Rndm()
+{
+   return (*this)();
+}
+
+double RanluxppCompatEngineStdRanlux24::operator()()
+{
+   return fImpl->NextRandomFloat();
+}
+
+uint64_t RanluxppCompatEngineStdRanlux24::IntRndm()
+{
+   return fImpl->NextRandomBits();
+}
+
+void RanluxppCompatEngineStdRanlux24::SetSeed(uint64_t seed)
+{
+   fImpl->SetSeedStd24(seed);
+}
+
+void RanluxppCompatEngineStdRanlux24::Skip(uint64_t n)
+{
+   fImpl->Skip(n);
+}
+
+
+RanluxppCompatEngineStdRanlux48::RanluxppCompatEngineStdRanlux48(uint64_t seed) : fImpl(new ImplType)
+{
+   this->SetSeed(seed);
+}
+
+RanluxppCompatEngineStdRanlux48::~RanluxppCompatEngineStdRanlux48() = default;
+
+double RanluxppCompatEngineStdRanlux48::Rndm()
+{
+   return (*this)();
+}
+
+double RanluxppCompatEngineStdRanlux48::operator()()
+{
+   return fImpl->NextRandomFloat();
+}
+
+uint64_t RanluxppCompatEngineStdRanlux48::IntRndm()
+{
+   return fImpl->NextRandomBits();
+}
+
+void RanluxppCompatEngineStdRanlux48::SetSeed(uint64_t seed)
+{
+   fImpl->SetSeedStd48(seed);
+}
+
+void RanluxppCompatEngineStdRanlux48::Skip(uint64_t n)
+{
+   fImpl->Skip(n);
+}
+
 } // end namespace Math
 } // end namespace ROOT
diff --git a/math/mathcore/test/RanluxppEngine.cxx b/math/mathcore/test/RanluxppEngine.cxx
index 258d90a055a7a..79c90237722e6 100644
--- a/math/mathcore/test/RanluxppEngine.cxx
+++ b/math/mathcore/test/RanluxppEngine.cxx
@@ -2,7 +2,7 @@
 // Author: Jonas Hahnfeld 11/2020
 
 /*************************************************************************
- * Copyright (C) 1995-2020, Rene Brun and Fons Rademakers.               *
+ * Copyright (C) 1995-2021, Rene Brun and Fons Rademakers.               *
  * All rights reserved.                                                  *
  *                                                                       *
  * For the licensing terms see $ROOTSYS/LICENSE.                         *
@@ -46,3 +46,460 @@ TEST(RanluxppEngine, random2048)
    EXPECT_EQ(rng.IntRndm(), 49145148745150);
    EXPECT_EQ(rng.Rndm(), 0.74670661284082484599);
 }
+
+TEST(RanluxppCompatEngineJames, P3)
+{
+   RanluxppCompatEngineJamesP3 rng(314159265);
+
+   // These values match: gsl_rng_ranlux, ranluxI_James, ranluxpp_James
+   EXPECT_EQ(rng.Rndm(), 0.53981816768646240234);
+   EXPECT_EQ(rng.Rndm(), 0.76155042648315429688);
+   EXPECT_EQ(rng.Rndm(), 0.06029939651489257812);
+   EXPECT_EQ(rng.Rndm(), 0.79600262641906738281);
+   EXPECT_EQ(rng.Rndm(), 0.30631220340728759766);
+
+   // Skip to the 24th number, right before the LCG is used to advance the state.
+   rng.Skip(18);
+   EXPECT_EQ(rng.Rndm(), 0.20569473505020141602);
+   // The LCG advances for the next call:
+   EXPECT_EQ(rng.Rndm(), 0.76727509498596191406);
+
+   // Skip to the 101st number.
+   rng.Skip(75);
+   EXPECT_EQ(rng.Rndm(), 0.43156743049621582031);
+   EXPECT_EQ(rng.Rndm(), 0.03774416446685791016);
+   EXPECT_EQ(rng.Rndm(), 0.24897110462188720703);
+   EXPECT_EQ(rng.Rndm(), 0.00147783756256103516);
+   EXPECT_EQ(rng.Rndm(), 0.90274453163146972656);
+}
+
+TEST(RanluxppCompatEngineJames, P4)
+{
+   RanluxppCompatEngineJamesP4 rng(314159265);
+
+   // These values match: gsl_rng_ranlux389
+   EXPECT_EQ(rng.Rndm(), 0.53981816768646240234);
+   EXPECT_EQ(rng.Rndm(), 0.76155042648315429688);
+   EXPECT_EQ(rng.Rndm(), 0.06029939651489257812);
+   EXPECT_EQ(rng.Rndm(), 0.79600262641906738281);
+   EXPECT_EQ(rng.Rndm(), 0.30631220340728759766);
+
+   // Skip to the 24th number, right before the LCG is used to advance the state.
+   rng.Skip(18);
+   EXPECT_EQ(rng.Rndm(), 0.20569473505020141602);
+   // The LCG advances for the next call:
+   EXPECT_EQ(rng.Rndm(), 0.84534603357315063477);
+
+   // Skip to the 101st number.
+   rng.Skip(75);
+   EXPECT_EQ(rng.Rndm(), 0.67576026916503906250);
+   EXPECT_EQ(rng.Rndm(), 0.90395343303680419922);
+   EXPECT_EQ(rng.Rndm(), 0.31414842605590820312);
+   EXPECT_EQ(rng.Rndm(), 0.98801732063293457031);
+   EXPECT_EQ(rng.Rndm(), 0.93221199512481689453);
+}
+
+TEST(RanluxppCompatEngineGslRanlxs, ranlxs0)
+{
+   RanluxppCompatEngineGslRanlxs0 rng(314159265);
+
+   // These values match: gsl_rng_ranlxs0
+   EXPECT_EQ(rng.Rndm(), 0.95476531982421875000);
+   EXPECT_EQ(rng.Rndm(), 0.10175001621246337891);
+   EXPECT_EQ(rng.Rndm(), 0.03923547267913818359);
+   EXPECT_EQ(rng.Rndm(), 0.23141473531723022461);
+   EXPECT_EQ(rng.Rndm(), 0.56545680761337280273);
+
+   // Skip to the 24th number, right before the LCG is used to advance the state.
+   rng.Skip(18);
+   EXPECT_EQ(rng.Rndm(), 0.66594201326370239258);
+   // The LCG advances for the next call:
+   EXPECT_EQ(rng.Rndm(), 0.08081126213073730469);
+
+   // Skip to the 101st number.
+   rng.Skip(75);
+   EXPECT_EQ(rng.Rndm(), 0.74328583478927612305);
+   EXPECT_EQ(rng.Rndm(), 0.79350239038467407227);
+   EXPECT_EQ(rng.Rndm(), 0.09384918212890625000);
+   EXPECT_EQ(rng.Rndm(), 0.00877797603607177734);
+   EXPECT_EQ(rng.Rndm(), 0.81286895275115966797);
+}
+
+TEST(RanluxppCompatEngineGslRanlxs, ranlxs0_default)
+{
+   RanluxppCompatEngineGslRanlxs0 rng;
+   EXPECT_EQ(rng.Rndm(), 0.32085895538330078125);
+
+   rng.SetSeed(0);
+   EXPECT_EQ(rng.Rndm(), 0.32085895538330078125);
+}
+
+TEST(RanluxppCompatEngineGslRanlxs, ranlxs1)
+{
+   RanluxppCompatEngineGslRanlxs1 rng(314159265);
+
+   // These values match: gsl_rng_ranlxs1
+   EXPECT_EQ(rng.Rndm(), 0.64368855953216552734);
+   EXPECT_EQ(rng.Rndm(), 0.67669576406478881836);
+   EXPECT_EQ(rng.Rndm(), 0.21001303195953369141);
+   EXPECT_EQ(rng.Rndm(), 0.27618223428726196289);
+   EXPECT_EQ(rng.Rndm(), 0.04699170589447021484);
+
+   // Skip to the 24th number, right before the LCG is used to advance the state.
+   rng.Skip(18);
+   EXPECT_EQ(rng.Rndm(), 0.81624960899353027344);
+   // The LCG advances for the next call:
+   EXPECT_EQ(rng.Rndm(), 0.52857619524002075195);
+
+   // Skip to the 101st number.
+   rng.Skip(75);
+   EXPECT_EQ(rng.Rndm(), 0.76676791906356811523);
+   EXPECT_EQ(rng.Rndm(), 0.28538894653320312500);
+   EXPECT_EQ(rng.Rndm(), 0.30258500576019287109);
+   EXPECT_EQ(rng.Rndm(), 0.38395434617996215820);
+   EXPECT_EQ(rng.Rndm(), 0.50305640697479248047);
+}
+
+TEST(RanluxppCompatEngineGslRanlxs, ranlxs1_default)
+{
+   RanluxppCompatEngineGslRanlxs1 rng;
+   EXPECT_EQ(rng.Rndm(), 0.06963491439819335938);
+
+   rng.SetSeed(0);
+   EXPECT_EQ(rng.Rndm(), 0.06963491439819335938);
+}
+
+TEST(RanluxppCompatEngineGslRanlxs, ranlxs2)
+{
+   RanluxppCompatEngineGslRanlxs2 rng(314159265);
+
+   // These values match: gsl_rng_ranlxs2
+   EXPECT_EQ(rng.Rndm(), 0.84824979305267333984);
+   EXPECT_EQ(rng.Rndm(), 0.75232243537902832031);
+   EXPECT_EQ(rng.Rndm(), 0.12018418312072753906);
+   EXPECT_EQ(rng.Rndm(), 0.05532860755920410156);
+   EXPECT_EQ(rng.Rndm(), 0.05795300006866455078);
+
+   // Skip to the 24th number, right before the LCG is used to advance the state.
+   rng.Skip(18);
+   EXPECT_EQ(rng.Rndm(), 0.13871461153030395508);
+   // The LCG advances for the next call:
+   EXPECT_EQ(rng.Rndm(), 0.97801673412322998047);
+
+   // Skip to the 101st number.
+   rng.Skip(75);
+   EXPECT_EQ(rng.Rndm(), 0.92746645212173461914);
+   EXPECT_EQ(rng.Rndm(), 0.82626664638519287109);
+   EXPECT_EQ(rng.Rndm(), 0.77763950824737548828);
+   EXPECT_EQ(rng.Rndm(), 0.49001514911651611328);
+   EXPECT_EQ(rng.Rndm(), 0.88770687580108642578);
+}
+
+TEST(RanluxppCompatEngineGslRanlxs, ranlxs2_default)
+{
+   RanluxppCompatEngineGslRanlxs2 rng;
+   EXPECT_EQ(rng.Rndm(), 0.53008824586868286133);
+
+   rng.SetSeed(0);
+   EXPECT_EQ(rng.Rndm(), 0.53008824586868286133);
+}
+
+TEST(RanluxppCompatEngineGslRanlxd, ranlxd1)
+{
+   RanluxppCompatEngineGslRanlxd1 rng(314159265);
+
+   // These values match: gsl_rng_ranlxd1
+   EXPECT_EQ(rng.Rndm(), 0.32330420393203596063);
+   EXPECT_EQ(rng.Rndm(), 0.72381776078560733367);
+   EXPECT_EQ(rng.Rndm(), 0.88817512439535306612);
+   EXPECT_EQ(rng.Rndm(), 0.04598644245910676887);
+   EXPECT_EQ(rng.Rndm(), 0.76110447268111158792);
+
+   // Skip to the 12th number, right before the LCG is used to advance the state.
+   rng.Skip(6);
+   EXPECT_EQ(rng.Rndm(), 0.18375035143688478456);
+   // The LCG advances for the next call:
+   EXPECT_EQ(rng.Rndm(), 0.83312931792853106572);
+
+   // Skip to the 101st number.
+   rng.Skip(87);
+   EXPECT_EQ(rng.Rndm(), 0.49329437415119770094);
+   EXPECT_EQ(rng.Rndm(), 0.37550852748539753634);
+   EXPECT_EQ(rng.Rndm(), 0.93543254436396239271);
+   EXPECT_EQ(rng.Rndm(), 0.43517686324045001811);
+   EXPECT_EQ(rng.Rndm(), 0.77751776122982363404);
+}
+
+TEST(RanluxppCompatEngineGslRanlxd, ranlxd1_default)
+{
+   RanluxppCompatEngineGslRanlxd1 rng;
+   EXPECT_EQ(rng.Rndm(), 0.83451879245814453157);
+
+   rng.SetSeed(0);
+   EXPECT_EQ(rng.Rndm(), 0.83451879245814453157);
+}
+
+TEST(RanluxppCompatEngineGslRanlxd, ranlxd2)
+{
+   RanluxppCompatEngineGslRanlxd2 rng(314159265);
+
+   // These values match: gsl_rng_ranlxd2
+   EXPECT_EQ(rng.Rndm(), 0.28994560412829528673);
+   EXPECT_EQ(rng.Rndm(), 0.06729010717305428102);
+   EXPECT_EQ(rng.Rndm(), 0.69404482039860582177);
+   EXPECT_EQ(rng.Rndm(), 0.56285566362747729841);
+   EXPECT_EQ(rng.Rndm(), 0.34655505440137091000);
+
+   // Skip to the 12th number, right before the LCG is used to advance the state.
+   rng.Skip(6);
+   EXPECT_EQ(rng.Rndm(), 0.12507187350916737500);
+   // The LCG advances for the next call:
+   EXPECT_EQ(rng.Rndm(), 0.97552452698574398937);
+
+   // Skip to the 101st number.
+   rng.Skip(87);
+   EXPECT_EQ(rng.Rndm(), 0.85118391727546338643);
+   EXPECT_EQ(rng.Rndm(), 0.94509863457767551154);
+   EXPECT_EQ(rng.Rndm(), 0.29384155017066149185);
+   EXPECT_EQ(rng.Rndm(), 0.48683495244473462549);
+   EXPECT_EQ(rng.Rndm(), 0.70125558306756730076);
+}
+
+TEST(RanluxppCompatEngineGslRanlxd, ranlxd2_default)
+{
+   RanluxppCompatEngineGslRanlxd2 rng;
+   EXPECT_EQ(rng.Rndm(), 0.07725383918716843823);
+
+   rng.SetSeed(0);
+   EXPECT_EQ(rng.Rndm(), 0.07725383918716843823);
+}
+
+TEST(RanluxppCompatEngineLuescherRanlxs, ranlxs0)
+{
+   RanluxppCompatEngineLuescherRanlxs0 rng(314159265);
+
+   // These values match: ranlxs0
+   EXPECT_EQ(rng.Rndm(), 0.65404480695724487305);
+   EXPECT_EQ(rng.Rndm(), 0.08226150274276733398);
+   EXPECT_EQ(rng.Rndm(), 0.29213166236877441406);
+   EXPECT_EQ(rng.Rndm(), 0.50535500049591064453);
+   EXPECT_EQ(rng.Rndm(), 0.95475935935974121094);
+
+   // Skip to the 92nd number, right before the LCG is used to advance the state.
+   rng.Skip(87);
+   EXPECT_EQ(rng.Rndm(), 0.63037145137786865234);
+   EXPECT_EQ(rng.Rndm(), 0.82800912857055664062);
+   EXPECT_EQ(rng.Rndm(), 0.31651723384857177734);
+   EXPECT_EQ(rng.Rndm(), 0.09057545661926269531);
+   // The LCG advances for the next calls:
+   EXPECT_EQ(rng.Rndm(), 0.06262010335922241211);
+   EXPECT_EQ(rng.Rndm(), 0.39436388015747070312);
+   EXPECT_EQ(rng.Rndm(), 0.61765056848526000977);
+   EXPECT_EQ(rng.Rndm(), 0.14016568660736083984);
+
+   // Skip to the 401st number.
+   rng.Skip(300);
+   EXPECT_EQ(rng.Rndm(), 0.21224951744079589844);
+   EXPECT_EQ(rng.Rndm(), 0.28637069463729858398);
+   EXPECT_EQ(rng.Rndm(), 0.01129972934722900391);
+   EXPECT_EQ(rng.Rndm(), 0.62236660718917846680);
+   EXPECT_EQ(rng.Rndm(), 0.11493819952011108398);
+}
+
+TEST(RanluxppCompatEngineLuescherRanlxs, ranlxs1)
+{
+   RanluxppCompatEngineLuescherRanlxs1 rng(314159265);
+
+   // These values match: ranlxs1
+   EXPECT_EQ(rng.Rndm(), 0.90099400281906127930);
+   EXPECT_EQ(rng.Rndm(), 0.76765865087509155273);
+   EXPECT_EQ(rng.Rndm(), 0.22530400753021240234);
+   EXPECT_EQ(rng.Rndm(), 0.83992105722427368164);
+   EXPECT_EQ(rng.Rndm(), 0.59816044569015502930);
+
+   // Skip to the 92nd number, right before the LCG is used to advance the state.
+   rng.Skip(87);
+   EXPECT_EQ(rng.Rndm(), 0.64669114351272583008);
+   EXPECT_EQ(rng.Rndm(), 0.46657902002334594727);
+   EXPECT_EQ(rng.Rndm(), 0.12610912322998046875);
+   EXPECT_EQ(rng.Rndm(), 0.10862994194030761719);
+   // The LCG advances for the next calls:
+   EXPECT_EQ(rng.Rndm(), 0.01416563987731933594);
+   EXPECT_EQ(rng.Rndm(), 0.58082711696624755859);
+   EXPECT_EQ(rng.Rndm(), 0.38216590881347656250);
+   EXPECT_EQ(rng.Rndm(), 0.91653412580490112305);
+
+   // Skip to the 401st number.
+   rng.Skip(300);
+   EXPECT_EQ(rng.Rndm(), 0.76539427042007446289);
+   EXPECT_EQ(rng.Rndm(), 0.11502689123153686523);
+   EXPECT_EQ(rng.Rndm(), 0.40491354465484619141);
+   EXPECT_EQ(rng.Rndm(), 0.96093446016311645508);
+   EXPECT_EQ(rng.Rndm(), 0.38819086551666259766);
+}
+
+TEST(RanluxppCompatEngineLuescherRanlxs, ranlxs2)
+{
+   RanluxppCompatEngineLuescherRanlxs2 rng(314159265);
+
+   // These values match: ranlxs2
+   EXPECT_EQ(rng.Rndm(), 0.60732543468475341797);
+   EXPECT_EQ(rng.Rndm(), 0.74212568998336791992);
+   EXPECT_EQ(rng.Rndm(), 0.76778668165206909180);
+   EXPECT_EQ(rng.Rndm(), 0.56459045410156250000);
+   EXPECT_EQ(rng.Rndm(), 0.51524215936660766602);
+
+   // Skip to the 92nd number, right before the LCG is used to advance the state.
+   rng.Skip(87);
+   EXPECT_EQ(rng.Rndm(), 0.07774782180786132812);
+   EXPECT_EQ(rng.Rndm(), 0.12600058317184448242);
+   EXPECT_EQ(rng.Rndm(), 0.56134593486785888672);
+   EXPECT_EQ(rng.Rndm(), 0.36321890354156494141);
+   // The LCG advances for the next calls:
+   EXPECT_EQ(rng.Rndm(), 0.23446798324584960938);
+   EXPECT_EQ(rng.Rndm(), 0.42847990989685058594);
+   EXPECT_EQ(rng.Rndm(), 0.21235740184783935547);
+   EXPECT_EQ(rng.Rndm(), 0.30497443675994873047);
+
+   // Skip to the 401st number.
+   rng.Skip(300);
+   EXPECT_EQ(rng.Rndm(), 0.17799735069274902344);
+   EXPECT_EQ(rng.Rndm(), 0.23861807584762573242);
+   EXPECT_EQ(rng.Rndm(), 0.65686619281768798828);
+   EXPECT_EQ(rng.Rndm(), 0.39222949743270874023);
+   EXPECT_EQ(rng.Rndm(), 0.45217937231063842773);
+}
+
+TEST(RanluxppCompatEngineLuescherRanlxd, ranlxd1)
+{
+   RanluxppCompatEngineLuescherRanlxd1 rng(314159265);
+
+   // These values match: ranlxd1
+   EXPECT_EQ(rng.Rndm(), 0.40183950697007020381);
+   EXPECT_EQ(rng.Rndm(), 0.75729181403551848462);
+   EXPECT_EQ(rng.Rndm(), 0.44293039330603889425);
+   EXPECT_EQ(rng.Rndm(), 0.27412463825007193918);
+   EXPECT_EQ(rng.Rndm(), 0.90543407471131232001);
+
+   // Skip to the 44th number, right before the LCG is used to advance the state.
+   rng.Skip(39);
+   EXPECT_EQ(rng.Rndm(), 0.35330884918042571030);
+   EXPECT_EQ(rng.Rndm(), 0.53342096246839787455);
+   EXPECT_EQ(rng.Rndm(), 0.87389084649229076263);
+   EXPECT_EQ(rng.Rndm(), 0.89137004735275837675);
+   // The LCG advances for the next calls:
+   EXPECT_EQ(rng.Rndm(), 0.92005646450699529737);
+   EXPECT_EQ(rng.Rndm(), 0.83437278691344118897);
+   EXPECT_EQ(rng.Rndm(), 0.96002403910338784954);
+   EXPECT_EQ(rng.Rndm(), 0.19537819235692666098);
+
+   // Skip to the 401st number.
+   rng.Skip(348);
+   EXPECT_EQ(rng.Rndm(), 0.74482704516052322674);
+   EXPECT_EQ(rng.Rndm(), 0.36930523933338221809);
+   EXPECT_EQ(rng.Rndm(), 0.36572707642276824913);
+   EXPECT_EQ(rng.Rndm(), 0.91548098968585378543);
+   EXPECT_EQ(rng.Rndm(), 0.55814623941527585771);
+}
+
+TEST(RanluxppCompatEngineLuescherRanlxd, ranlxd2)
+{
+   RanluxppCompatEngineLuescherRanlxd2 rng(314159265);
+
+   // These values match: ranlxd2
+   EXPECT_EQ(rng.Rndm(), 0.52702589450092673928);
+   EXPECT_EQ(rng.Rndm(), 0.14545363991469173470);
+   EXPECT_EQ(rng.Rndm(), 0.79180414135301191436);
+   EXPECT_EQ(rng.Rndm(), 0.07733853533090595533);
+   EXPECT_EQ(rng.Rndm(), 0.14180497711929263005);
+
+   // Skip to the 44th number, right before the LCG is used to advance the state.
+   rng.Skip(39);
+   EXPECT_EQ(rng.Rndm(), 0.18603867668678475411);
+   EXPECT_EQ(rng.Rndm(), 0.13778587858302770996);
+   EXPECT_EQ(rng.Rndm(), 0.70244055389935411426);
+   EXPECT_EQ(rng.Rndm(), 0.90056757149901045523);
+   // The LCG advances for the next calls:
+   EXPECT_EQ(rng.Rndm(), 0.93044152534523050235);
+   EXPECT_EQ(rng.Rndm(), 0.99330367160287025285);
+   EXPECT_EQ(rng.Rndm(), 0.63528838564139178402);
+   EXPECT_EQ(rng.Rndm(), 0.61276861283659656010);
+
+   // Skip to the 401st number.
+   rng.Skip(348);
+   EXPECT_EQ(rng.Rndm(), 0.37774257837078550892);
+   EXPECT_EQ(rng.Rndm(), 0.20167640311612444748);
+   EXPECT_EQ(rng.Rndm(), 0.34101141228811471251);
+   EXPECT_EQ(rng.Rndm(), 0.35947589472738883387);
+   EXPECT_EQ(rng.Rndm(), 0.51160837322844798791);
+}
+
+TEST(RanluxppCompatEngineStdRanlux24, compare)
+{
+   RanluxppCompatEngineStdRanlux24 rng(314159265);
+
+   // These values match: std::ranlux24
+   EXPECT_EQ(rng.IntRndm(), 6389521);
+   EXPECT_EQ(rng.IntRndm(), 1245860);
+   EXPECT_EQ(rng.IntRndm(), 9047089);
+   EXPECT_EQ(rng.IntRndm(), 5613314);
+   EXPECT_EQ(rng.IntRndm(), 15388463);
+
+   // Skip to the 23rd number, right before the LCG is used to advance the state.
+   rng.Skip(17);
+   EXPECT_EQ(rng.IntRndm(), 14135596);
+   // The LCG advances for the next call:
+   EXPECT_EQ(rng.IntRndm(), 1842057);
+
+   // Skip to the 101st number.
+   rng.Skip(76);
+   EXPECT_EQ(rng.IntRndm(), 7894321);
+   EXPECT_EQ(rng.IntRndm(), 2634015);
+   EXPECT_EQ(rng.IntRndm(), 12134196);
+   EXPECT_EQ(rng.IntRndm(), 15231589);
+   EXPECT_EQ(rng.IntRndm(), 11032869);
+}
+
+TEST(RanluxppCompatEngineStdRanlux24, default)
+{
+   RanluxppCompatEngineStdRanlux24 rng;
+
+   // Skip to the 10000th number, which is specified by the standard.
+   rng.Skip(9999);
+   EXPECT_EQ(rng.IntRndm(), 9901578);
+}
+
+TEST(RanluxppCompatEngineStdRanlux48, compare)
+{
+   RanluxppCompatEngineStdRanlux48 rng(314159265);
+
+   // These values match: std::ranlux48
+   EXPECT_EQ(rng.IntRndm(), 2902733192977);
+   EXPECT_EQ(rng.IntRndm(), 183625379875889);
+   EXPECT_EQ(rng.IntRndm(), 164401649471280);
+   EXPECT_EQ(rng.IntRndm(), 158572192479531);
+   EXPECT_EQ(rng.IntRndm(), 227024878504710);
+
+   // Skip to the 11th number, right before the LCG is used to advance the state.
+   rng.Skip(5);
+   EXPECT_EQ(rng.IntRndm(), 137631957902972);
+   // The LCG advances for the next call:
+   EXPECT_EQ(rng.IntRndm(), 122438241205867);
+
+   // Skip to the 101st number.
+   rng.Skip(88);
+   EXPECT_EQ(rng.IntRndm(), 155713325118081);
+   EXPECT_EQ(rng.IntRndm(), 75203262258561);
+   EXPECT_EQ(rng.IntRndm(), 164155826303104);
+   EXPECT_EQ(rng.IntRndm(), 58159697115827);
+   EXPECT_EQ(rng.IntRndm(), 89006261856016);
+}
+
+TEST(RanluxppCompatEngineStdRanlux48, default)
+{
+   RanluxppCompatEngineStdRanlux48 rng;
+
+   // Skip to the 10000th number, which is specified by the standard.
+   rng.Skip(9999);
+   EXPECT_EQ(rng.IntRndm(), 249142670248501);
+}