From c6a6499a82fbbabac8eaebf3107332c110484e19 Mon Sep 17 00:00:00 2001 From: Jonas Hahnfeld Date: Wed, 9 Jun 2021 12:32:12 +0200 Subject: [PATCH] RANLUX++: Add compatibility engines MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit These engines can be used to obtain the same sequences of numbers as RANLUX generators using recursive subtract-with-borrow steps, but with enhanced performance. Apart from the choice of parameters, the main difference between the various implementations is the way of seeding the initial state of the generator. This commit includes engines for compatibility with: * the original implementation by Fred James, with parameters for - luxury level 3 (p = 223), also matching gsl_rng_ranlux - luxury level 4 (p = 389), also matching gsl_rng_ranlux389 producing floating point numbers from 24 bits of randomness; * the family of generators using a second-generation version of the RANLUX algorithm as implemented in the GNU Scientific Library: - gsl_rng_ranlxs[012] using 24 bits per floating point number, and - gsl_rng_ranlxd[12] using 48 bits per floating point number; * the implementation by Martin Lüscher written in C that uses four states per generator; similar to GSL, there are ranlxs[012] with 24 bits per number and ranlxd[12] with 48 bits per number; and * the generators std::ranlux{24,48} defined by the C++ standard. The values in the tests were extracted directly from the mentioned implementations, showing that the LCG implementation is equivalent to the RANLUX algorithm. I am not adding compatibility engines for CLHEP because its semantics are very weird: While CLHEP::RanluxEngine::setSeed yields the same sequences as the original implementation by James, the seed is treated differently when passed as an argument to the constructor. --- math/mathcore/inc/Math/RanluxppEngine.h | 247 ++++++++- math/mathcore/src/RanluxppEngineImpl.cxx | 609 ++++++++++++++++++++++- math/mathcore/test/RanluxppEngine.cxx | 459 ++++++++++++++++- 3 files changed, 1300 insertions(+), 15 deletions(-) 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 +template class RanluxppEngineImpl; template class RanluxppEngine final : public TRandomEngine { private: - std::unique_ptr> fImpl; + using ImplType = RanluxppEngineImpl<48, p>; + std::unique_ptr 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 +class RanluxppCompatEngineJames final : public TRandomEngine { + +private: + using ImplType = RanluxppEngineImpl<24, p>; + std::unique_ptr 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 +class RanluxppCompatEngineGslRanlxs final : public TRandomEngine { + +private: + using ImplType = RanluxppEngineImpl<24, p>; + std::unique_ptr 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 +class RanluxppCompatEngineGslRanlxd final : public TRandomEngine { + +private: + using ImplType = RanluxppEngineImpl<48, p>; + std::unique_ptr 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 +class RanluxppCompatEngineLuescherImpl; + +/// Compatibility engine for Lüscher's ranlxs implementation written in C. +template +class RanluxppCompatEngineLuescherRanlxs final : public TRandomEngine { + +private: + using ImplType = RanluxppCompatEngineLuescherImpl<24, p>; + std::unique_ptr 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 +class RanluxppCompatEngineLuescherRanlxd final : public TRandomEngine { + +private: + using ImplType = RanluxppCompatEngineLuescherImpl<48, p>; + std::unique_ptr 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 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 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, , m) +// >>> print_hex(kA) + template 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 +template class RanluxppEngineImpl { + // Needs direct access to private members to initialize its four states. + friend class RanluxppCompatEngineLuescherImpl; 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

::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 -RanluxppEngine

::RanluxppEngine(uint64_t seed) : fImpl(new RanluxppEngineImpl<48, p>) +RanluxppEngine

::RanluxppEngine(uint64_t seed) : fImpl(new ImplType) { - fImpl->SetSeed(seed); + this->SetSeed(seed); } template @@ -212,7 +428,7 @@ uint64_t RanluxppEngine

::IntRndm() template void RanluxppEngine

::SetSeed(uint64_t seed) { - fImpl->SetSeed(seed); + fImpl->SetSeedSibidanov(seed); } template @@ -224,5 +440,376 @@ void RanluxppEngine

::Skip(uint64_t n) template class RanluxppEngine<24>; template class RanluxppEngine<2048>; + +template +RanluxppCompatEngineJames

::RanluxppCompatEngineJames(uint64_t seed) : fImpl(new ImplType) +{ + this->SetSeed(seed); +} + +template +RanluxppCompatEngineJames

::~RanluxppCompatEngineJames() = default; + +template +double RanluxppCompatEngineJames

::Rndm() +{ + return (*this)(); +} + +template +double RanluxppCompatEngineJames

::operator()() +{ + return fImpl->NextRandomFloat(); +} + +template +uint64_t RanluxppCompatEngineJames

::IntRndm() +{ + return fImpl->NextRandomBits(); +} + +template +void RanluxppCompatEngineJames

::SetSeed(uint64_t seed) +{ + fImpl->SetSeedJames(seed); +} + +template +void RanluxppCompatEngineJames

::Skip(uint64_t n) +{ + fImpl->Skip(n); +} + +template class RanluxppCompatEngineJames<223>; +template class RanluxppCompatEngineJames<389>; + + +template +RanluxppCompatEngineGslRanlxs

::RanluxppCompatEngineGslRanlxs(uint64_t seed) : fImpl(new ImplType) +{ + this->SetSeed(seed); +} + +template +RanluxppCompatEngineGslRanlxs

::~RanluxppCompatEngineGslRanlxs() = default; + +template +double RanluxppCompatEngineGslRanlxs

::Rndm() +{ + return (*this)(); +} + +template +double RanluxppCompatEngineGslRanlxs

::operator()() +{ + return fImpl->NextRandomFloat(); +} + +template +uint64_t RanluxppCompatEngineGslRanlxs

::IntRndm() +{ + return fImpl->NextRandomBits(); +} + +template +void RanluxppCompatEngineGslRanlxs

::SetSeed(uint64_t seed) +{ + fImpl->SetSeedGsl(seed, /*ranlxd=*/false); +} + +template +void RanluxppCompatEngineGslRanlxs

::Skip(uint64_t n) +{ + fImpl->Skip(n); +} + +template class RanluxppCompatEngineGslRanlxs<218>; +template class RanluxppCompatEngineGslRanlxs<404>; +template class RanluxppCompatEngineGslRanlxs<794>; + + +template +RanluxppCompatEngineGslRanlxd

::RanluxppCompatEngineGslRanlxd(uint64_t seed) : fImpl(new ImplType) +{ + this->SetSeed(seed); +} + +template +RanluxppCompatEngineGslRanlxd

::~RanluxppCompatEngineGslRanlxd() = default; + +template +double RanluxppCompatEngineGslRanlxd

::Rndm() +{ + return (*this)(); +} + +template +double RanluxppCompatEngineGslRanlxd

::operator()() +{ + return fImpl->NextRandomFloat(); +} + +template +uint64_t RanluxppCompatEngineGslRanlxd

::IntRndm() +{ + return fImpl->NextRandomBits(); +} + +template +void RanluxppCompatEngineGslRanlxd

::SetSeed(uint64_t seed) +{ + fImpl->SetSeedGsl(seed, /*ranlxd=*/true); +} + +template +void RanluxppCompatEngineGslRanlxd

::Skip(uint64_t n) +{ + fImpl->Skip(n); +} + +template class RanluxppCompatEngineGslRanlxd<404>; +template class RanluxppCompatEngineGslRanlxd<794>; + + +template +class RanluxppCompatEngineLuescherImpl { + +private: + RanluxppEngineImpl 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 +RanluxppCompatEngineLuescherRanlxs

::RanluxppCompatEngineLuescherRanlxs(uint64_t seed) : fImpl(new ImplType) +{ + this->SetSeed(seed); +} + +template +RanluxppCompatEngineLuescherRanlxs

::~RanluxppCompatEngineLuescherRanlxs() = default; + +template +double RanluxppCompatEngineLuescherRanlxs

::Rndm() +{ + return (*this)(); +} + +template +double RanluxppCompatEngineLuescherRanlxs

::operator()() +{ + return fImpl->NextRandomFloat(); +} + +template +uint64_t RanluxppCompatEngineLuescherRanlxs

::IntRndm() +{ + return fImpl->NextRandomBits(); +} + +template +void RanluxppCompatEngineLuescherRanlxs

::SetSeed(uint64_t seed) +{ + fImpl->SetSeed(seed, /*ranlxd=*/false); +} + +template +void RanluxppCompatEngineLuescherRanlxs

::Skip(uint64_t n) +{ + fImpl->Skip(n); +} + +template class RanluxppCompatEngineLuescherRanlxs<218>; +template class RanluxppCompatEngineLuescherRanlxs<404>; +template class RanluxppCompatEngineLuescherRanlxs<794>; + + +template +RanluxppCompatEngineLuescherRanlxd

::RanluxppCompatEngineLuescherRanlxd(uint64_t seed) : fImpl(new ImplType) +{ + this->SetSeed(seed); +} + +template +RanluxppCompatEngineLuescherRanlxd

::~RanluxppCompatEngineLuescherRanlxd() = default; + +template +double RanluxppCompatEngineLuescherRanlxd

::Rndm() +{ + return (*this)(); +} + +template +double RanluxppCompatEngineLuescherRanlxd

::operator()() +{ + return fImpl->NextRandomFloat(); +} + +template +uint64_t RanluxppCompatEngineLuescherRanlxd

::IntRndm() +{ + return fImpl->NextRandomBits(); +} + +template +void RanluxppCompatEngineLuescherRanlxd

::SetSeed(uint64_t seed) +{ + fImpl->SetSeed(seed, /*ranlxd=*/true); +} + +template +void RanluxppCompatEngineLuescherRanlxd

::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); +}