RANLUX++: Add compatibility engines #8383
Conversation
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.
|
Starting build on |
|
To elaborate a bit on CLHEP: CLHEP::RanluxEngine r;
r.setSeed(314159265);but directly passing the seed to the constructor à la diff --git a/math/mathcore/src/RanluxppEngineImpl.cxx b/math/mathcore/src/RanluxppEngineImpl.cxx
index 100f8d8638..bbf508a6a8 100644
--- a/math/mathcore/src/RanluxppEngineImpl.cxx
+++ b/math/mathcore/src/RanluxppEngineImpl.cxx
@@ -219,13 +219,14 @@ public:
// 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;
+ int64_t seed = s & 0xffffff;
auto next = [&]() {
const int a = 0xd1a4, b = 0x9c4e, c = 0x2fb3, d = 0x7fffffab;
+ int64_t oldSeed = seed;
int64_t k = seed / a;
seed = b * (seed - k * a) - k * c ;
if (seed < 0) seed += d;
- return seed & 0xffffff;
+ return oldSeed & 0xffffff;
};
// Iteration is reversed because the first number from the MCG goes to theThat would add compatibility for the constructor, but leave no way to call |
|
Can you motivate why we should include those in ROOT's interface? I understand the motivation for testing! I'm sure you have a good reason to also expose them, I'd just like to see the reasons :-) |
|
@Axel-Naumann yes, testing is one of the motivations, in particular continuous testing to prevent future regressions (now we can check against an external implementation, instead of just copying the current values and declaring them "known-good"). The other reason, and why I think this might provide benefit for users, is performance: The original RANLUX implementation by James (at least its implementation in GSL) needs 40 seconds to sum 1 million numbers at luxury level 3, Now we could argue that all users should switch to |
|
ping @lmoneta |
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.
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:
p = 223), also matchinggsl_rng_ranluxp = 389), also matchinggsl_rng_ranlux389producing floating point numbers from 24 bits of randomness;gsl_rng_ranlxs[012]using 24 bits per floating point number, andgsl_rng_ranlxd[12]using 48 bits per floating point number;ranlxs[012]with 24 bits per number andranlxd[12]with 48 bits per number; andstd::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::setSeedyields the same sequences as the original implementation by James, the seed is treated differently when passed as an argument to the constructor.