|
RandomLib
1.10
|
|
RandomLib
1.10
|
This page contains the following sections
Here are simple steps to convert a code using rand() to Random:
These examples use Random::Global which is a globally defined Random object. In many cases this suffices. However, if different parts of your code need independent random number streams, you can create your own Random objects. In these cases, you might want to vary the seeds used to initialize the separate streams in a systematic way and, in that case, you can seed Random with a vector. When you pass Random objects to other routines you should usually (always?) pass them by reference, Random&, to ensure that that the change in the state of the generator is seen by the parent routine.
Do not pass a Random object to std::generate to fill an vector with random integers. Even though operator()() is defined to produce an random integer in [0,232), generate makes a copy of its operator argument. Thus
results in a and b having the same contents and r.Count() == 0. You should instead pass the Random object by reference as in
Alternatively (and more flexibly) you can define function objects which accept the Random by reference in the constructor as follows:
Now, you can use
C++11 offers a better way of doing this using lambda expressions. These allow you to write function objects "in line" (rather than as some small once-used class in another part of the file. See RandomLambda.cpp.
On the other hand, std::random_shuffle does pass its operator argument by reference. Thus a vector can be shuffled with
This shuffles because operator()(unsigned n) is defined to produce an random integer in [0,n). This shuffling can result in all possible permutations of vectors of lengths up to 2000 (because 2000! < 219937), whereas the built-in random number generator (invoked when the last argument to random_shuffle is omitted) can typically only produce all the permutations of 12 or fewer objects. This method of shuffling only works if the number of elements being shuffled is less than 2^32. On 64-bit computers, this condition can be violated. In that case, use SRandom64 as the basic random number class, or else pass a function object which can accept a unsigned long long argument; for an example see RandomLambda.cpp.
When saving the state of a program in a restart file, it is usually necessary to save only the Random objects with Save(...). (In fact, Random contains no additional state beyond the state of RandomGenerator. And thus Save just calls the underlying RandomGenerator::Save.)
RandomSelect has state which is derivable from its input weights. However, neither it nor the other distributions contain state which depends on the Random argument to operator()(Random& r) (which is a const member function).
For speed and for better control over round-off errors, the real routines FixedX are preferred over FloatX. In a few cases, FloatX gives better results, e.g., in the implementation of ExponentialDistribution where it provides finer granularity in the results.
Fixed, FixedU, FixedN, FixedW, FixedS are all obtained by rounding an ideal uniform deviate and so can all be used to sample periodic intervals uniformly. If possible, avoid using FixedC and FixedO since they can introduce bias into your simulations. Instead of FixedC, consider FixedN. Instead of FixedO, consider FixedU (to avoid 0) or shift interval to (−1/2,1/2) and use FixedS. You can also generate results in (0,1) by invoking FixedS with a smaller precision and shifting the result, for example
This library defines a STATIC_ASSERT to check template parameters and the values of some constants at compile time. For example, this will prevent you from requesting 100 bits of accuracy in a float, e.g., Fixed<float,100>(). The error message you get from the compiler may not be very informative. However if you look at the corresponding line of source code in the header file, you should be able to figure out the problem.
Avoid using calling RandomLib twice in one expression, e.g.,
because you don't know which random number will of computed first. This means that changing compilers or even just changing the optimization level might lead to large changes in your results even if the random seed is the same. In this case, you could instead use
Similarly replace
by
RandomLib provides 8 typedefs offering you easy access to 4 different random generators:
Thus {M,S}RandomGenerator{32,64} gives access to the 32-bit and 64-bit versions of MT19937 and SFMT19937 generators.
Normally, two additional typedefs are provided
However if the preprocessor symbol RANDOMLIB_DEFAULT_GENERATOR is defined to be one of {M,S}RandomGenerator{32,64} when compiling files that include Random.hpp, then RandomGenerator is defined to be this instead. This can be defined by supplying -DRANDOMLIB_DEFAULT_GENERATOR=MRandomGenerator32 on the command line for the compiler or by specifying the value when invoking make, e.g., make RANDOMLIB_DEFAULT_GENERATOR=MRandomGenerator32 RandomExample.
In normal use, all the functions of RandomLib should be accessed via the Random class (or occasionally via the underlying generator RandomGenerator).
The SIMD-oriented Fast Mersenne Twister random number generator, SFMT19937 was developed in 2006 as an improvement on the MT19937. By and large these two generators have very similar properties and can be regarded as strong enough for nearly all applications. However the SFMT19937 generator does have some advantages:
In addition, SFMT19937 adopted an improved scheme for converting the seed into the state (via MixerSFMT) and, by default, this is used for the MT19937 generator in this library. Because of these advantages, Random is typedef'ed to the SFMT19937 generator, by default.
The 32-bit and 64-bit versions are comparable in strength. Indeed with the SFMT19937 generator these are essentially the same (the underlying algorithm manipulates 128-bit words in both cases). For example, we have
The implementations of these is portable across 32-bit and 64-bit architectures. Thus the choice of between these is probably best made on the basis of the speed on the target platform. A glance at the timing data in the next section shows that the 64-bit versions are about the same speed on a 64-bit machine (x86_64). However on a 64-bit machine, random routines which can consume 64-bits random results in one piece, e.g., r.Integer<unsigned long long>() and r.Fixed<double>, are faster with 64-bit versions of the generator. In practice, the 32-bit versions are more likely to perform well on a wide range of CPUs. However, it's likely that the 64-bit version will be a better choice in a few years (particularly if double-precision floating point random numbers are used).
I recommend against mixing generators within a single code. This will merely result in a more complicated code. However, while you are experimenting with different generators, you should print out the type of generator being used, given by Name(), together with the seed. For example:
Quick summary: stick with the default generator SRandom32 accessed through Random. In certain cases, switching to the 64-bit version SRandom64 might be advisable.
The following times were obtained on a Linux system running Fedora 12, kernel version 2.6.32, and compiling with g++ version 4.4.4 with optimization flags "-O3 -funroll-loops", on an Intel system. Here r is a Random object, n is a unsigned variable with value 0 (but the compiler doesn't know its value), N is a large positive integer, and all times are in ns (unless another time unit is given)
| operation | Intel Xeon, x86_64 2.66GHz (SSE2) | |
|---|---|---|
| SRandom32 | SRandom64 | |
| std::rand() | 10 | |
| r() | 2.0 | 2.7 |
| r.Integer<unsigned>() | 2.3 | 2.6 |
| r.Integer<unsigned long long>() | 4.3 | 2.7 |
| r.Integer<unsigned,6>() | 2.0 | 2.9 |
| r.Integer<unsigned>(52u) | 3.0 | 3.7 |
| r.Integer<unsigned>(52u+n) | 10 | 11 |
| r.Fixed<float>() | 3.1 | 3.5 |
| r.Fixed<double>() | 5.6 | 3.4 |
| r.Float<float>() | 17 | 18 |
| r.Float<double>() | 18 | 17 |
| r.Prob<float>(0.28f) | 9.1 | 9.3 |
| r.Prob<double>(0.28) | 8.0 | 8.6 |
| NormalDistribution<float>()(r) | 35 | 35 |
| NormalDistribution<double>()(r) | 44 | 43 |
| SeedWord() | 1.1 ms | |
| SeedVector() | 19 us | |
| r.Reset(), r.SetCount(0) | 6.9 us | |
| r.StepCount(N) | 0.63 N | 1.3 N |
| r.StepCount(-N) | 1.1 N | 2.3 N |
These timing figures were produced by RandomTime.cpp. You should run this yourself to determine the times relevant to your system. It is easy to adapt RandomTime.cpp to produce timings for MRandom{32,64}; however on SSE2 systems, these will be noticeably slower than the SFMT19937 generators.
Loop unrolling is critical in the performance of, e.g., Fixed<double>(). If your compiler doesn't unroll the loops in Fixed<RealType,p>(), you can provide specializations with the loops unrolled by hand. Much of the speed of this implementation comes from extensive use of inlined procedures. This also make the timing results sensitive to the context in which Random is called. Compare the results for r.Integer<unsigned>(52u) (typical if drawing a card from a deck) and r.Integer<unsigned>(52u+n) (as typical called by std::random_shuffle). In the first case the compiler can precompute some of the variables used resulting is a much faster execution. The best way of gauging the speed is to time or to profile your own application.
The time quoted for r.Reset(0), r.SetCount() gives the time to reset the generator and to convert the seed to the initial state. Note that r.Reseed(...) merely stores the seed and does not update the state and r.Reset() also similarly does not reinitialize the generator state. The random generator state is produced from the seed by calling the mixer when the first random number is requested or when SetCount is called.
Beware of the compiler optimizing too much code away when doing timing studies. Timing this section of code
will usually result is an unrealistically short time (by up to a factor of three!) because the compiler skips over a lot of the computation (the tempering of the random results and all the real multiplications). You can prevent the compiler from "cheating" in this way by using the results from the timed functions. Thus
More simply, you can store the results in a variable which is declared volatile
The SelfTest() method of the random generators checks the that the correct results are obtained stepping the generator both forwards and backwards. This throws an exception if the an error is detected. E.g.,
RandomExample.cpp listed in Getting started includes tests of all the generators.
In addition, RandomCoverage.cpp includes code to check the various generators using their authors' test cases. This checks the current implementations of MT19937 and SFMT19937 against their original implementations. (The original 64-bit implementation of MT19937 used a vector of 64-bit integers for the seed. So for this check above it is necessary to recast these as a vector of 32-bit integers,by inserting 0 elements.)
1.8.6
