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
RandomLib::Random r; std::vector<unsigned> a(10); std::vector<unsigned> b(10); std::generate(a.begin(), a.end(), r); std::generate(b.begin(), b.end(), r);
results in a and b having the same contents and r.Count() == 0. You should instead pass the Random object by reference as in
std::generate<std::vector<unsigned>::iterator, RandomLib::Random&> (a.begin(), a.end(), r); std::generate<std::vector<unsigned>::iterator, RandomLib::Random&> (b.begin(), b.end(), r);
Alternatively (and more flexibly) you can define function objects which accept the Random by reference in the constructor as follows:
template<typename IntType = int> class RandomInteger { private: RandomLib::Random& _r; const IntType _m, _n; public: RandomInteger(RandomLib::Random& r, IntType m, IntType n) : _r(r), _m(m), _n(n) {} IntType operator()() { return _r.IntegerC<IntType>(_m, _n); } }; template<typename RealType = double> class RandomNormal { private: RandomLib::Random& _r; const RandomLib::NormalDistribution<RealType> _n; const RealType _mean, _sigma; public: RandomNormal(RandomLib::Random& r, RealType mean = RealType(0), RealType sigma = RealType(1)) : _r(r), _n(RandomLib::NormalDistribution<RealType>()) , _mean(mean), _sigma(sigma) {} RealType operator()() { return _n(_r, _mean, _sigma); } };
Now, you can use
RandomLib::Random r; std::vector<int> a(10); // Fill with integers in [-10,10] std::generate(a.begin(), a.end(), RandomInteger<>(r,-10,10)); std::vector<double> b(10); // Fill with normal deviates std::generate(b.begin(), b.end(), RandomNormal<>(r,1.0,2.0));
On the other hand, std::random_shuffle does pass its operator argument by reference. Thus a vector can be shuffled with
RandomLib::Random r; std::vector<unsigned> a(100); for (unsigned i = 0; i < 100; ++i) a[i] = i; // create a random permutation of [0,100) std::random_shuffle(a.begin(), a.end(), r);
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.
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
RandomLib::Random r; double y = // result in (0,1) r.FixedS<double, std::numeric_limits<double>::digits - 1>() + 0.5;
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.
RandomLib provides 4 typedefs:
Thus {M,S}RandomGenerator{32,64} gives access to the 32-bit and 64-bit versions of MT19937 and SFMT19937 generators.
Normally, an additional typedefs is provided
However if the preprocessor symbol 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 -DDEFAULT_GENERATOR=MRandomGenerator32 on the command line for the compiler or by specifying the value when invoking make, e.g., make DEFAULT_GENERATOR=MRandomGenerator32 RandomExample.
Finally, we define
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, after
we have a.Ran64() == b.Ran64().
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 slower than the 32-bit version on a 32-bit architecture (i686) and is 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 practise, 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 withing a single code. This will merely result in a more complicated code. However, while you are experimenting with different generators, I suggest printing out the type of generator being used, given by Name(), together with the seed. For example:
RandomLib::RandomCanonical<SRandomGenerator64> r; std::cout << "Using " << r.Name() << " with seed " << r.SeedString() << std::endl;
The following times were obtained on Linux systems running Fedora 7 or 8, kernel version 2.6.23, and compiling with g++ version 4.1.2 with optimization flags "-O3 -funroll-loops -finline-functions -fomit-frame-pointer". The random objects {M,S}Random{32,64} are here abbreviated {MT,SFMT}{32,64}. 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 | Pentium III, i686, 864MHz | Intel Xeon, x86_64 3.6GHz (SSE2) | ||||||
|---|---|---|---|---|---|---|---|---|
| MT32 | MT64 | SFMT32 | SFMT64 | MT32 | MT64 | SFMT32 | SFMT64 | |
| std::rand() | 86 | 18 | ||||||
| r() | 25 | 66 | 19 | 48 | 6.0 | 8.8 | 2.5 | 3.8 |
| r.Integer<unsigned>() | 25 | 64 | 20 | 47 | 6.0 | 8.8 | 2.5 | 3.8 |
| r.Integer<unsigned long long>() | 51 | 66 | 43 | 48 | 12 | 8.8 | 5.0 | 3.8 |
| r.Integer<unsigned,6>() | 26 | 65 | 20 | 48 | 6.1 | 9.3 | 2.7 | 3.9 |
| r.Integer<unsigned>(52u) | 35 | 73 | 27 | 54 | 8.2 | 11 | 4.4 | 6.1 |
| r.Integer<unsigned>(52u+n) | 98 | 140 | 110 | 140 | 29 | 31 | 29 | 31 |
| r.Fixed<float>() | 47 | 92 | 46 | 80 | 9.6 | 12 | 6.0 | 7.5 |
| r.Fixed<double>() | 110 | 90 | 110 | 80 | 24 | 13 | 19 | 7.9 |
| r.Float<float>() | 100 | 160 | 90 | 140 | 37 | 44 | 32 | 35 |
| r.Float<double>() | 170 | 170 | 150 | 140 | 52 | 43 | 45 | 34 |
| r.Prob<float>(0.28f) | 67 | 110 | 62 | 94 | 18 | 22 | 15 | 16 |
| r.Prob<double>(0.28) | 74 | 120 | 71 | 100 | 19 | 22 | 16 | 17 |
| NormalDistribution<float>()(r) | 250 | 370 | 240 | 340 | 80 | 90 | 73 | 80 |
| NormalDistribution<double>()(r) | 380 | 370 | 370 | 340 | 110 | 87 | 92 | 80 |
| SeedWord() | 6.1 ms | 1.4 ms | ||||||
| SeedVector() | 80 us | 25 us | ||||||
| r.Reset(), r.SetCount(0) | 45 us | 10 us | ||||||
| r.StepCount(N) | 8.8 N | 34 N | 13 N | 40 N | 2.3 N | 3.1 N | 1.3 N | 2.6 N |
| r.StepCount(-N) | 13 N | 40 N | 15 N | 43 N | 3.1 N | 4.7 N | 1.9 N | 3.9 N |
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(), 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
RandomLib::Random r; const size_t n = 100000000; double d = 0; for (size_t i = n; i; --i) d += r.Fixed<double>(); std::cout << "Sum: " << d << std::endl;
More simply, you can store the results in a variable which is declared volatile
RandomLib::Random r; const size_t n = 100000000; volatile double d; for (size_t i = n; i; --i) d = r.Fixed<double>();
As of 2006-11-03, at least one implementor of the Mersenne Twister makes a bogus speed claim for his implementations because of a poorly designed timing test:
There's nothing the matter with this implementation. It's just that it isn't as fast as claimed.
Compile the following code (and Random.cpp) with -DRANDOM_LEGACY
using namespace RandomLib; #if RANDOM_LEGACY { RandomEngine<MT19937<Random_u32>,MixerMT0<Random_u32> > s("0x123,0x234,0x345,0x456"); s.SetCount(999); uint32_t x = s(); std::cout << "Generator: " << s.Name() << "\n" << "Seed: " << s.SeedString() << "\n" << s.Count() << "th result: " << x << " (should be 3460025646)"<< "\n\n"; } { RandomEngine<MT19937<Random_u64>,MixerMT0<Random_u64> > s("0x12345,0,0x23456,0,0x34567,0,0x45678,0"); s.SetCount(999); uint64_t x = s(); std::cout << "Generator: " << s.Name() << "\n" << "Seed: " << s.SeedString() << "\n" << s.Count() << "th result: " << x << " (should be 994412663058993407)"<< "\n\n"; } #endif { SRandomGenerator32 s("0x1234,0x5678,0x9abc,0xdef0"); s.SetCount(999); uint32_t x = s(); std::cout << "Generator: " << s.Name() << "\n" << "Seed: " << s.SeedString() << "\n" << s.Count() << "th result: " << x << " (should be 788493625)"<< "\n\n"; } { SRandomGenerator64 s("5,4,3,2,1"); s.SetCount(999); uint64_t x = s(); std::cout << "Generator: " << s.Name() << "\n" << "Seed: " << s.SeedString() << "\n" << s.Count() << "th result: " << x << " (should be 13356980519185762498)"<< "\n\n"; }
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 the check above it is necessary to recast these as a vector of 32-bit integers (by inserting 0 elements). This should produce the following output:
Generator: RandomEngine<MT19937<Random_u32>,MixerMT0<Random_u32>> Seed: [291,564,837,1110] 1000th result: 3460025646 (should be 3460025646) Generator: RandomEngine<MT19937<Random_u64>,MixerMT0<Random_u64>> Seed: [74565,0,144470,0,214375,0,284280,0] 1000th result: 994412663058993407 (should be 994412663058993407) Generator: RandomEngine<SFMT19937<Random_u32>,MixerSFMT> Seed: [4660,22136,39612,57072] 1000th result: 788493625 (should be 788493625) Generator: RandomEngine<SFMT19937<Random_u64>,MixerSFMT> Seed: [5,4,3,2,1] 1000th result: 13356980519185762498 (should be 13356980519185762498)
1.7.5