Chris@16: /* boost random/mersenne_twister.hpp header file Chris@16: * Chris@16: * Copyright Jens Maurer 2000-2001 Chris@16: * Copyright Steven Watanabe 2010 Chris@16: * Distributed under the Boost Software License, Version 1.0. (See Chris@16: * accompanying file LICENSE_1_0.txt or copy at Chris@16: * http://www.boost.org/LICENSE_1_0.txt) Chris@16: * Chris@16: * See http://www.boost.org for most recent version including documentation. Chris@16: * Chris@16: * $Id: mersenne_twister.hpp 74867 2011-10-09 23:13:31Z steven_watanabe $ Chris@16: * Chris@16: * Revision history Chris@16: * 2001-02-18 moved to individual header files Chris@16: */ Chris@16: Chris@16: #ifndef BOOST_RANDOM_MERSENNE_TWISTER_HPP Chris@16: #define BOOST_RANDOM_MERSENNE_TWISTER_HPP Chris@16: Chris@16: #include Chris@16: #include Chris@16: #include Chris@16: #include Chris@16: #include Chris@16: #include Chris@16: #include Chris@16: #include Chris@16: #include Chris@16: #include Chris@16: #include Chris@16: Chris@16: namespace boost { Chris@16: namespace random { Chris@16: Chris@16: /** Chris@16: * Instantiations of class template mersenne_twister_engine model a Chris@16: * \pseudo_random_number_generator. It uses the algorithm described in Chris@16: * Chris@16: * @blockquote Chris@16: * "Mersenne Twister: A 623-dimensionally equidistributed uniform Chris@16: * pseudo-random number generator", Makoto Matsumoto and Takuji Nishimura, Chris@16: * ACM Transactions on Modeling and Computer Simulation: Special Issue on Chris@16: * Uniform Random Number Generation, Vol. 8, No. 1, January 1998, pp. 3-30. Chris@16: * @endblockquote Chris@16: * Chris@16: * @xmlnote Chris@16: * The boost variant has been implemented from scratch and does not Chris@16: * derive from or use mt19937.c provided on the above WWW site. However, it Chris@16: * was verified that both produce identical output. Chris@16: * @endxmlnote Chris@16: * Chris@16: * The seeding from an integer was changed in April 2005 to address a Chris@16: * weakness. Chris@16: * Chris@16: * The quality of the generator crucially depends on the choice of the Chris@16: * parameters. User code should employ one of the sensibly parameterized Chris@16: * generators such as \mt19937 instead. Chris@16: * Chris@16: * The generator requires considerable amounts of memory for the storage of Chris@16: * its state array. For example, \mt11213b requires about 1408 bytes and Chris@16: * \mt19937 requires about 2496 bytes. Chris@16: */ Chris@16: template Chris@16: class mersenne_twister_engine Chris@16: { Chris@16: public: Chris@16: typedef UIntType result_type; Chris@16: BOOST_STATIC_CONSTANT(std::size_t, word_size = w); Chris@16: BOOST_STATIC_CONSTANT(std::size_t, state_size = n); Chris@16: BOOST_STATIC_CONSTANT(std::size_t, shift_size = m); Chris@16: BOOST_STATIC_CONSTANT(std::size_t, mask_bits = r); Chris@16: BOOST_STATIC_CONSTANT(UIntType, xor_mask = a); Chris@16: BOOST_STATIC_CONSTANT(std::size_t, tempering_u = u); Chris@16: BOOST_STATIC_CONSTANT(UIntType, tempering_d = d); Chris@16: BOOST_STATIC_CONSTANT(std::size_t, tempering_s = s); Chris@16: BOOST_STATIC_CONSTANT(UIntType, tempering_b = b); Chris@16: BOOST_STATIC_CONSTANT(std::size_t, tempering_t = t); Chris@16: BOOST_STATIC_CONSTANT(UIntType, tempering_c = c); Chris@16: BOOST_STATIC_CONSTANT(std::size_t, tempering_l = l); Chris@16: BOOST_STATIC_CONSTANT(UIntType, initialization_multiplier = f); Chris@16: BOOST_STATIC_CONSTANT(UIntType, default_seed = 5489u); Chris@16: Chris@16: // backwards compatibility Chris@16: BOOST_STATIC_CONSTANT(UIntType, parameter_a = a); Chris@16: BOOST_STATIC_CONSTANT(std::size_t, output_u = u); Chris@16: BOOST_STATIC_CONSTANT(std::size_t, output_s = s); Chris@16: BOOST_STATIC_CONSTANT(UIntType, output_b = b); Chris@16: BOOST_STATIC_CONSTANT(std::size_t, output_t = t); Chris@16: BOOST_STATIC_CONSTANT(UIntType, output_c = c); Chris@16: BOOST_STATIC_CONSTANT(std::size_t, output_l = l); Chris@16: Chris@16: // old Boost.Random concept requirements Chris@16: BOOST_STATIC_CONSTANT(bool, has_fixed_range = false); Chris@16: Chris@16: Chris@16: /** Chris@16: * Constructs a @c mersenne_twister_engine and calls @c seed(). Chris@16: */ Chris@16: mersenne_twister_engine() { seed(); } Chris@16: Chris@16: /** Chris@16: * Constructs a @c mersenne_twister_engine and calls @c seed(value). Chris@16: */ Chris@16: BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(mersenne_twister_engine, Chris@16: UIntType, value) Chris@16: { seed(value); } Chris@16: template mersenne_twister_engine(It& first, It last) Chris@16: { seed(first,last); } Chris@16: Chris@16: /** Chris@16: * Constructs a mersenne_twister_engine and calls @c seed(gen). Chris@16: * Chris@16: * @xmlnote Chris@16: * The copy constructor will always be preferred over Chris@16: * the templated constructor. Chris@16: * @endxmlnote Chris@16: */ Chris@16: BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(mersenne_twister_engine, Chris@16: SeedSeq, seq) Chris@16: { seed(seq); } Chris@16: Chris@16: // compiler-generated copy ctor and assignment operator are fine Chris@16: Chris@16: /** Calls @c seed(default_seed). */ Chris@16: void seed() { seed(default_seed); } Chris@16: Chris@16: /** Chris@16: * Sets the state x(0) to v mod 2w. Then, iteratively, Chris@16: * sets x(i) to Chris@16: * (i + f * (x(i-1) xor (x(i-1) rshift w-2))) mod 2w Chris@16: * for i = 1 .. n-1. x(n) is the first value to be returned by operator(). Chris@16: */ Chris@16: BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(mersenne_twister_engine, UIntType, value) Chris@16: { Chris@16: // New seeding algorithm from Chris@16: // http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html Chris@16: // In the previous versions, MSBs of the seed affected only MSBs of the Chris@16: // state x[]. Chris@16: const UIntType mask = (max)(); Chris@16: x[0] = value & mask; Chris@16: for (i = 1; i < n; i++) { Chris@16: // See Knuth "The Art of Computer Programming" Chris@16: // Vol. 2, 3rd ed., page 106 Chris@16: x[i] = (f * (x[i-1] ^ (x[i-1] >> (w-2))) + i) & mask; Chris@16: } Chris@16: } Chris@16: Chris@16: /** Chris@16: * Seeds a mersenne_twister_engine using values produced by seq.generate(). Chris@16: */ Chris@16: BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(mersenne_twister_engine, SeeqSeq, seq) Chris@16: { Chris@16: detail::seed_array_int(seq, x); Chris@16: i = n; Chris@16: Chris@16: // fix up the state if it's all zeroes. Chris@16: if((x[0] & (~static_cast(0) << r)) == 0) { Chris@16: for(std::size_t j = 1; j < n; ++j) { Chris@16: if(x[j] != 0) return; Chris@16: } Chris@16: x[0] = static_cast(1) << (w-1); Chris@16: } Chris@16: } Chris@16: Chris@16: /** Sets the state of the generator using values from an iterator range. */ Chris@16: template Chris@16: void seed(It& first, It last) Chris@16: { Chris@16: detail::fill_array_int(first, last, x); Chris@16: i = n; Chris@16: Chris@16: // fix up the state if it's all zeroes. Chris@16: if((x[0] & (~static_cast(0) << r)) == 0) { Chris@16: for(std::size_t j = 1; j < n; ++j) { Chris@16: if(x[j] != 0) return; Chris@16: } Chris@16: x[0] = static_cast(1) << (w-1); Chris@16: } Chris@16: } Chris@16: Chris@16: /** Returns the smallest value that the generator can produce. */ Chris@16: static result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () Chris@16: { return 0; } Chris@16: /** Returns the largest value that the generator can produce. */ Chris@16: static result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () Chris@16: { return boost::low_bits_mask_t::sig_bits; } Chris@16: Chris@16: /** Produces the next value of the generator. */ Chris@16: result_type operator()(); Chris@16: Chris@16: /** Fills a range with random values */ Chris@16: template Chris@16: void generate(Iter first, Iter last) Chris@16: { detail::generate_from_int(*this, first, last); } Chris@16: Chris@16: /** Chris@16: * Advances the state of the generator by @c z steps. Equivalent to Chris@16: * Chris@16: * @code Chris@16: * for(unsigned long long i = 0; i < z; ++i) { Chris@16: * gen(); Chris@16: * } Chris@16: * @endcode Chris@16: */ Chris@16: void discard(boost::uintmax_t z) Chris@16: { Chris@16: for(boost::uintmax_t j = 0; j < z; ++j) { Chris@16: (*this)(); Chris@16: } Chris@16: } Chris@16: Chris@16: #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS Chris@16: /** Writes a mersenne_twister_engine to a @c std::ostream */ Chris@16: template Chris@16: friend std::basic_ostream& Chris@16: operator<<(std::basic_ostream& os, Chris@16: const mersenne_twister_engine& mt) Chris@16: { Chris@16: mt.print(os); Chris@16: return os; Chris@16: } Chris@16: Chris@16: /** Reads a mersenne_twister_engine from a @c std::istream */ Chris@16: template Chris@16: friend std::basic_istream& Chris@16: operator>>(std::basic_istream& is, Chris@16: mersenne_twister_engine& mt) Chris@16: { Chris@16: for(std::size_t j = 0; j < mt.state_size; ++j) Chris@16: is >> mt.x[j] >> std::ws; Chris@16: // MSVC (up to 7.1) and Borland (up to 5.64) don't handle the template Chris@16: // value parameter "n" available from the class template scope, so use Chris@16: // the static constant with the same value Chris@16: mt.i = mt.state_size; Chris@16: return is; Chris@16: } Chris@16: #endif Chris@16: Chris@16: /** Chris@16: * Returns true if the two generators are in the same state, Chris@16: * and will thus produce identical sequences. Chris@16: */ Chris@16: friend bool operator==(const mersenne_twister_engine& x, Chris@16: const mersenne_twister_engine& y) Chris@16: { Chris@16: if(x.i < y.i) return x.equal_imp(y); Chris@16: else return y.equal_imp(x); Chris@16: } Chris@16: Chris@16: /** Chris@16: * Returns true if the two generators are in different states. Chris@16: */ Chris@16: friend bool operator!=(const mersenne_twister_engine& x, Chris@16: const mersenne_twister_engine& y) Chris@16: { return !(x == y); } Chris@16: Chris@16: private: Chris@16: /// \cond show_private Chris@16: Chris@16: void twist(); Chris@16: Chris@16: /** Chris@16: * Does the work of operator==. This is in a member function Chris@16: * for portability. Some compilers, such as msvc 7.1 and Chris@16: * Sun CC 5.10 can't access template parameters or static Chris@16: * members of the class from inline friend functions. Chris@16: * Chris@16: * requires i <= other.i Chris@16: */ Chris@16: bool equal_imp(const mersenne_twister_engine& other) const Chris@16: { Chris@16: UIntType back[n]; Chris@16: std::size_t offset = other.i - i; Chris@16: for(std::size_t j = 0; j + offset < n; ++j) Chris@16: if(x[j] != other.x[j+offset]) Chris@16: return false; Chris@16: rewind(&back[n-1], offset); Chris@16: for(std::size_t j = 0; j < offset; ++j) Chris@16: if(back[j + n - offset] != other.x[j]) Chris@16: return false; Chris@16: return true; Chris@16: } Chris@16: Chris@16: /** Chris@16: * Does the work of operator<<. This is in a member function Chris@16: * for portability. Chris@16: */ Chris@16: template Chris@16: void print(std::basic_ostream& os) const Chris@16: { Chris@16: UIntType data[n]; Chris@16: for(std::size_t j = 0; j < i; ++j) { Chris@16: data[j + n - i] = x[j]; Chris@16: } Chris@16: if(i != n) { Chris@16: rewind(&data[n - i - 1], n - i); Chris@16: } Chris@16: os << data[0]; Chris@16: for(std::size_t j = 1; j < n; ++j) { Chris@16: os << ' ' << data[j]; Chris@16: } Chris@16: } Chris@16: Chris@16: /** Chris@16: * Copies z elements of the state preceding x[0] into Chris@16: * the array whose last element is last. Chris@16: */ Chris@16: void rewind(UIntType* last, std::size_t z) const Chris@16: { Chris@16: const UIntType upper_mask = (~static_cast(0)) << r; Chris@16: const UIntType lower_mask = ~upper_mask; Chris@16: UIntType y0 = x[m-1] ^ x[n-1]; Chris@16: if(y0 & (static_cast(1) << (w-1))) { Chris@16: y0 = ((y0 ^ a) << 1) | 1; Chris@16: } else { Chris@16: y0 = y0 << 1; Chris@16: } Chris@16: for(std::size_t sz = 0; sz < z; ++sz) { Chris@16: UIntType y1 = Chris@16: rewind_find(last, sz, m-1) ^ rewind_find(last, sz, n-1); Chris@16: if(y1 & (static_cast(1) << (w-1))) { Chris@16: y1 = ((y1 ^ a) << 1) | 1; Chris@16: } else { Chris@16: y1 = y1 << 1; Chris@16: } Chris@16: *(last - sz) = (y0 & upper_mask) | (y1 & lower_mask); Chris@16: y0 = y1; Chris@16: } Chris@16: } Chris@16: Chris@16: /** Chris@16: * Given a pointer to the last element of the rewind array, Chris@16: * and the current size of the rewind array, finds an element Chris@16: * relative to the next available slot in the rewind array. Chris@16: */ Chris@16: UIntType Chris@16: rewind_find(UIntType* last, std::size_t size, std::size_t j) const Chris@16: { Chris@16: std::size_t index = (j + n - size + n - 1) % n; Chris@16: if(index < n - size) { Chris@16: return x[index]; Chris@16: } else { Chris@16: return *(last - (n - 1 - index)); Chris@16: } Chris@16: } Chris@16: Chris@16: /// \endcond Chris@16: Chris@16: // state representation: next output is o(x(i)) Chris@16: // x[0] ... x[k] x[k+1] ... x[n-1] represents Chris@16: // x(i-k) ... x(i) x(i+1) ... x(i-k+n-1) Chris@16: Chris@16: UIntType x[n]; Chris@16: std::size_t i; Chris@16: }; Chris@16: Chris@16: /// \cond show_private Chris@16: Chris@16: #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION Chris@16: // A definition is required even for integral static constants Chris@16: #define BOOST_RANDOM_MT_DEFINE_CONSTANT(type, name) \ Chris@16: template \ Chris@16: const type mersenne_twister_engine::name Chris@16: BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, word_size); Chris@16: BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, state_size); Chris@16: BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, shift_size); Chris@16: BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, mask_bits); Chris@16: BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, xor_mask); Chris@16: BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, tempering_u); Chris@16: BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, tempering_d); Chris@16: BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, tempering_s); Chris@16: BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, tempering_b); Chris@16: BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, tempering_t); Chris@16: BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, tempering_c); Chris@16: BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, tempering_l); Chris@16: BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, initialization_multiplier); Chris@16: BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, default_seed); Chris@16: BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, parameter_a); Chris@16: BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, output_u ); Chris@16: BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, output_s); Chris@16: BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, output_b); Chris@16: BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, output_t); Chris@16: BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, output_c); Chris@16: BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, output_l); Chris@16: BOOST_RANDOM_MT_DEFINE_CONSTANT(bool, has_fixed_range); Chris@16: #undef BOOST_RANDOM_MT_DEFINE_CONSTANT Chris@16: #endif Chris@16: Chris@16: template Chris@16: void Chris@16: mersenne_twister_engine::twist() Chris@16: { Chris@16: const UIntType upper_mask = (~static_cast(0)) << r; Chris@16: const UIntType lower_mask = ~upper_mask; Chris@16: Chris@16: const std::size_t unroll_factor = 6; Chris@16: const std::size_t unroll_extra1 = (n-m) % unroll_factor; Chris@16: const std::size_t unroll_extra2 = (m-1) % unroll_factor; Chris@16: Chris@16: // split loop to avoid costly modulo operations Chris@16: { // extra scope for MSVC brokenness w.r.t. for scope Chris@16: for(std::size_t j = 0; j < n-m-unroll_extra1; j++) { Chris@16: UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask); Chris@16: x[j] = x[j+m] ^ (y >> 1) ^ ((x[j+1]&1) * a); Chris@16: } Chris@16: } Chris@16: { Chris@16: for(std::size_t j = n-m-unroll_extra1; j < n-m; j++) { Chris@16: UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask); Chris@16: x[j] = x[j+m] ^ (y >> 1) ^ ((x[j+1]&1) * a); Chris@16: } Chris@16: } Chris@16: { Chris@16: for(std::size_t j = n-m; j < n-1-unroll_extra2; j++) { Chris@16: UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask); Chris@16: x[j] = x[j-(n-m)] ^ (y >> 1) ^ ((x[j+1]&1) * a); Chris@16: } Chris@16: } Chris@16: { Chris@16: for(std::size_t j = n-1-unroll_extra2; j < n-1; j++) { Chris@16: UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask); Chris@16: x[j] = x[j-(n-m)] ^ (y >> 1) ^ ((x[j+1]&1) * a); Chris@16: } Chris@16: } Chris@16: // last iteration Chris@16: UIntType y = (x[n-1] & upper_mask) | (x[0] & lower_mask); Chris@16: x[n-1] = x[m-1] ^ (y >> 1) ^ ((x[0]&1) * a); Chris@16: i = 0; Chris@16: } Chris@16: /// \endcond Chris@16: Chris@16: template Chris@16: inline typename Chris@16: mersenne_twister_engine::result_type Chris@16: mersenne_twister_engine::operator()() Chris@16: { Chris@16: if(i == n) Chris@16: twist(); Chris@16: // Step 4 Chris@16: UIntType z = x[i]; Chris@16: ++i; Chris@16: z ^= ((z >> u) & d); Chris@16: z ^= ((z << s) & b); Chris@16: z ^= ((z << t) & c); Chris@16: z ^= (z >> l); Chris@16: return z; Chris@16: } Chris@16: Chris@16: /** Chris@16: * The specializations \mt11213b and \mt19937 are from Chris@16: * Chris@16: * @blockquote Chris@16: * "Mersenne Twister: A 623-dimensionally equidistributed Chris@16: * uniform pseudo-random number generator", Makoto Matsumoto Chris@16: * and Takuji Nishimura, ACM Transactions on Modeling and Chris@16: * Computer Simulation: Special Issue on Uniform Random Number Chris@16: * Generation, Vol. 8, No. 1, January 1998, pp. 3-30. Chris@16: * @endblockquote Chris@16: */ Chris@16: typedef mersenne_twister_engine mt11213b; Chris@16: Chris@16: /** Chris@16: * The specializations \mt11213b and \mt19937 are from Chris@16: * Chris@16: * @blockquote Chris@16: * "Mersenne Twister: A 623-dimensionally equidistributed Chris@16: * uniform pseudo-random number generator", Makoto Matsumoto Chris@16: * and Takuji Nishimura, ACM Transactions on Modeling and Chris@16: * Computer Simulation: Special Issue on Uniform Random Number Chris@16: * Generation, Vol. 8, No. 1, January 1998, pp. 3-30. Chris@16: * @endblockquote Chris@16: */ Chris@16: typedef mersenne_twister_engine mt19937; Chris@16: Chris@16: #if !defined(BOOST_NO_INT64_T) && !defined(BOOST_NO_INTEGRAL_INT64_T) Chris@16: typedef mersenne_twister_engine mt19937_64; Chris@16: #endif Chris@16: Chris@16: /// \cond show_deprecated Chris@16: Chris@16: template Chris@16: class mersenne_twister : Chris@16: public mersenne_twister_engine Chris@16: { Chris@16: typedef mersenne_twister_engine base_type; Chris@16: public: Chris@16: mersenne_twister() {} Chris@16: BOOST_RANDOM_DETAIL_GENERATOR_CONSTRUCTOR(mersenne_twister, Gen, gen) Chris@16: { seed(gen); } Chris@16: BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(mersenne_twister, UIntType, val) Chris@16: { seed(val); } Chris@16: template Chris@16: mersenne_twister(It& first, It last) : base_type(first, last) {} Chris@16: void seed() { base_type::seed(); } Chris@16: BOOST_RANDOM_DETAIL_GENERATOR_SEED(mersenne_twister, Gen, gen) Chris@16: { Chris@16: detail::generator_seed_seq seq(gen); Chris@16: base_type::seed(seq); Chris@16: } Chris@16: BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(mersenne_twister, UIntType, val) Chris@16: { base_type::seed(val); } Chris@16: template Chris@16: void seed(It& first, It last) { base_type::seed(first, last); } Chris@16: }; Chris@16: Chris@16: /// \endcond Chris@16: Chris@16: } // namespace random Chris@16: Chris@16: using random::mt11213b; Chris@16: using random::mt19937; Chris@16: using random::mt19937_64; Chris@16: Chris@16: } // namespace boost Chris@16: Chris@16: BOOST_RANDOM_PTR_HELPER_SPEC(boost::mt11213b) Chris@16: BOOST_RANDOM_PTR_HELPER_SPEC(boost::mt19937) Chris@16: BOOST_RANDOM_PTR_HELPER_SPEC(boost::mt19937_64) Chris@16: Chris@16: #endif // BOOST_RANDOM_MERSENNE_TWISTER_HPP