diff --git a/include/boost/random/detail/seed.hpp b/include/boost/random/detail/seed.hpp index c090b18..55b2fa6 100644 --- a/include/boost/random/detail/seed.hpp +++ b/include/boost/random/detail/seed.hpp @@ -43,12 +43,19 @@ struct disable_constructor {}; template \ void seed(Generator& gen, typename ::boost::random::detail::disable_seed::type* = 0) +#define BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(Self, SeedSeq, seq) \ + template \ + explicit Self(SeedSeq& seq, typename ::boost::random::detail::disable_constructor::type* = 0) + +#define BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(Self, SeedSeq, seq) \ + template \ + void seed(SeedSeq& seq, typename ::boost::random::detail::disable_seed::type* = 0) + #define BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(Self, T, x) \ explicit Self(const T& x) #define BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(Self, T, x) \ void seed(const T& x) - } } } @@ -76,6 +83,24 @@ struct disable_constructor {}; template\ void boost_random_seed_impl(Generator& gen, ::boost::mpl::false_) +#define BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(Self, SeedSeq, seq) \ + Self(Self& other) { *this = other; } \ + Self(const Self& other) { *this = other; } \ + template \ + explicit Self(SeedSeq& seq) { \ + boost_random_constructor_impl(seq, ::boost::is_arithmetic());\ + } \ + template \ + void boost_random_constructor_impl(SeedSeq& seq, ::boost::mpl::false_) + +#define BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(Self, SeedSeq, seq) \ + template \ + void seed(SeedSeq& seq) { \ + boost_random_seed_impl(seq, ::boost::is_arithmetic()); \ + } \ + template \ + void boost_random_seed_impl(SeedSeq& seq, ::boost::mpl::false_) + #define BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(Self, T, x) \ explicit Self(const T& x) { boost_random_constructor_impl(x, ::boost::mpl::true_()); }\ void boost_random_constructor_impl(const T& x, ::boost::mpl::true_) diff --git a/include/boost/random/linear_congruential.hpp b/include/boost/random/linear_congruential.hpp index 01f12f6..ab9b2a3 100644 --- a/include/boost/random/linear_congruential.hpp +++ b/include/boost/random/linear_congruential.hpp @@ -160,6 +160,15 @@ public: _x = const_mod::mult_add(a, _x, c); return _x; } + + /** Fills a range with random values */ + template + void generate(Iter first, Iter last) + { + for(; first != last; ++first) { + *first = (*this)(); + } + } static bool validation(IntType x) { return val == x; } diff --git a/include/boost/random/mersenne_twister.hpp b/include/boost/random/mersenne_twister.hpp index 4215353..0f4e81a 100644 --- a/include/boost/random/mersenne_twister.hpp +++ b/include/boost/random/mersenne_twister.hpp @@ -1,6 +1,7 @@ /* boost random/mersenne_twister.hpp header file * * Copyright Jens Maurer 2000-2001 + * Copyright Steven Watanabe 2010 * Distributed under the Boost Software License, Version 1.0. (See * accompanying file LICENSE_1_0.txt or copy at * http://www.boost.org/LICENSE_1_0.txt) @@ -16,16 +17,10 @@ #ifndef BOOST_RANDOM_MERSENNE_TWISTER_HPP #define BOOST_RANDOM_MERSENNE_TWISTER_HPP -#include -#include // std::copy +#include #include #include -#include -#include -#include #include -#include -#include #include #include #include @@ -34,7 +29,7 @@ namespace boost { namespace random { /** - * Instantiations of class template mersenne_twister model a + * Instantiations of class template mersenne_twister_engine model a * \pseudo_random_number_generator. It uses the algorithm described in * * @blockquote @@ -61,307 +56,375 @@ namespace random { * its state array. For example, \mt11213b requires about 1408 bytes and * \mt19937 requires about 2496 bytes. */ -template -class mersenne_twister +template +class mersenne_twister_engine { public: - typedef UIntType result_type; - BOOST_STATIC_CONSTANT(int, word_size = w); - BOOST_STATIC_CONSTANT(int, state_size = n); - BOOST_STATIC_CONSTANT(int, shift_size = m); - BOOST_STATIC_CONSTANT(int, mask_bits = r); - BOOST_STATIC_CONSTANT(UIntType, parameter_a = a); - BOOST_STATIC_CONSTANT(int, output_u = u); - BOOST_STATIC_CONSTANT(int, output_s = s); - BOOST_STATIC_CONSTANT(UIntType, output_b = b); - BOOST_STATIC_CONSTANT(int, output_t = t); - BOOST_STATIC_CONSTANT(UIntType, output_c = c); - BOOST_STATIC_CONSTANT(int, output_l = l); - - BOOST_STATIC_CONSTANT(bool, has_fixed_range = false); + typedef UIntType result_type; + BOOST_STATIC_CONSTANT(std::size_t, word_size = w); + BOOST_STATIC_CONSTANT(std::size_t, state_size = n); + BOOST_STATIC_CONSTANT(std::size_t, shift_size = m); + BOOST_STATIC_CONSTANT(std::size_t, mask_bits = r); + BOOST_STATIC_CONSTANT(UIntType, xor_mask = a); + BOOST_STATIC_CONSTANT(std::size_t, tempering_u = u); + BOOST_STATIC_CONSTANT(UIntType, tempering_d = d); + BOOST_STATIC_CONSTANT(std::size_t, tempering_s = s); + BOOST_STATIC_CONSTANT(UIntType, tempering_b = b); + BOOST_STATIC_CONSTANT(std::size_t, tempering_t = t); + BOOST_STATIC_CONSTANT(UIntType, tempering_c = c); + BOOST_STATIC_CONSTANT(std::size_t, tempering_l = l); + BOOST_STATIC_CONSTANT(UIntType, initialization_multiplier = f); + BOOST_STATIC_CONSTANT(UIntType, default_seed = 5489u); - /** - * Constructs a @c mersenne_twister and calls @c seed(). - */ - mersenne_twister() { seed(); } + // backwards compatibility + BOOST_STATIC_CONSTANT(UIntType, parameter_a = a); + BOOST_STATIC_CONSTANT(std::size_t, output_u = u); + BOOST_STATIC_CONSTANT(std::size_t, output_s = s); + BOOST_STATIC_CONSTANT(UIntType, output_b = b); + BOOST_STATIC_CONSTANT(std::size_t, output_t = t); + BOOST_STATIC_CONSTANT(UIntType, output_c = c); + BOOST_STATIC_CONSTANT(std::size_t, output_l = l); + + // old Boost.Random concept requirements + BOOST_STATIC_CONSTANT(bool, has_fixed_range = false); - /** - * Constructs a @c mersenne_twister and calls @c seed(value). - */ - BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(mersenne_twister, UIntType, value) - { seed(value); } - template mersenne_twister(It& first, It last) { seed(first,last); } - /** - * Constructs a mersenne_twister and calls @c seed(gen). - * - * @xmlnote - * The copy constructor will always be preferred over - * the templated constructor. - * @endxmlnote - */ - BOOST_RANDOM_DETAIL_GENERATOR_CONSTRUCTOR(mersenne_twister, Generator, gen) - { seed(gen); } + /** + * Constructs a @c mersenne_twister_engine and calls @c seed(). + */ + mersenne_twister_engine() { seed(); } - // compiler-generated copy ctor and assignment operator are fine + /** + * Constructs a @c mersenne_twister_engine and calls @c seed(value). + */ + BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(mersenne_twister_engine, UIntType, value) + { seed(value); } + template mersenne_twister_engine(It& first, It last) + { seed(first,last); } - /** Calls @c seed(result_type(5489)). */ - void seed() { seed(UIntType(5489)); } + /** + * Constructs a mersenne_twister_engine and calls @c seed(gen). + * + * @xmlnote + * The copy constructor will always be preferred over + * the templated constructor. + * @endxmlnote + */ + BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(mersenne_twister_engine, SeedSeq, seq) + { seed(seq); } - /** - * Sets the state x(0) to v mod 2w. Then, iteratively, - * sets x(i) to (i + 1812433253 * (x(i-1) xor (x(i-1) rshift w-2))) mod 2w - * for i = 1 .. n-1. x(n) is the first value to be returned by operator(). - */ - BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(mersenne_twister, UIntType, value) - { - // New seeding algorithm from - // http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html - // In the previous versions, MSBs of the seed affected only MSBs of the - // state x[]. - const UIntType mask = ~0u; - x[0] = value & mask; - for (i = 1; i < n; i++) { - // See Knuth "The Art of Computer Programming" Vol. 2, 3rd ed., page 106 - x[i] = (1812433253UL * (x[i-1] ^ (x[i-1] >> (w-2))) + i) & mask; + // compiler-generated copy ctor and assignment operator are fine + + /** Calls @c seed(default_seed). */ + void seed() { seed(default_seed); } + + /** + * Sets the state x(0) to v mod 2w. Then, iteratively, + * sets x(i) to + * (i + f * (x(i-1) xor (x(i-1) rshift w-2))) mod 2w + * for i = 1 .. n-1. x(n) is the first value to be returned by operator(). + */ + BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(mersenne_twister_engine, UIntType, value) + { + // New seeding algorithm from + // http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html + // In the previous versions, MSBs of the seed affected only MSBs of the + // state x[]. + const UIntType mask = (max)(); + x[0] = value & mask; + for (i = 1; i < n; i++) { + // See Knuth "The Art of Computer Programming" Vol. 2, 3rd ed., page 106 + x[i] = (f * (x[i-1] ^ (x[i-1] >> (w-2))) + i) & mask; + } } - } + + /** + * Seeds a mersenne_twister_engine using values produced by seq.generate(). + */ + BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(mersenne_twister_engine, SeeqSeq, seq) + { + boost::uint_least32_t storage[((w+31)/32) * n]; + seq.generate(&storage[0], &storage[0] + ((w+31)/32) * n); + for(std::size_t j = 0; j < n; j++) { + UIntType val = 0; + for(std::size_t k = 0; k < (w+31)/32; ++k) { + val += storage[(w+31)/32*j + k] << 32*k; + } + x[j] = val; + } + i = n; - /** - * Sets the state of this mersenne_twister to the values - * returned by n invocations of gen. - * - * Complexity: Exactly n invocations of gen. - */ - BOOST_RANDOM_DETAIL_GENERATOR_SEED(mersenne_twister, Generator, gen) - { -#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS - BOOST_STATIC_ASSERT(!std::numeric_limits::is_signed); -#endif - // I could have used std::generate_n, but it takes "gen" by value - for(int j = 0; j < n; j++) - x[j] = gen(); - i = n; - } + // fix up the state if it's all zeroes. + if((x[0] & (~static_cast(0) << r)) == 0) { + for(std::size_t j = 1; i < n; ++j) { + if(x[j] != 0) return; + } + x[0] = static_cast(1) << (w-1); + } + } - template - void seed(It& first, It last) - { - int j; - for(j = 0; j < n && first != last; ++j, ++first) - x[j] = *first; - i = n; - if(first == last && j < n) - throw std::invalid_argument("mersenne_twister::seed"); - } + /** Sets the state of the generator using values from an iterator range. */ + template + void seed(It& first, It last) + { + std::size_t j; + for(j = 0; j < n && first != last; ++j, ++first) + x[j] = *first; + i = n; + if(first == last && j < n) + throw std::invalid_argument("mersenne_twister_engine::seed"); + } - result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; } - result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const - { - // avoid "left shift count >= with of type" warning - result_type res = 0; - for(int j = 0; j < w; ++j) - res |= (1u << j); - return res; - } + /** Returns the smallest value that the generator can produce. */ + static result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () + { return 0; } + /** Returns the largest value that the generator can produce. */ + static result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () + { + // avoid "left shift count >= width of type" warning + result_type res = 0; + for(std::size_t j = 0; j < w; ++j) + res |= (static_cast(1) << j); + return res; + } + + /** Produces the next value of the generator. */ + result_type operator()(); - result_type operator()(); - static bool validation(result_type v) { return val == v; } + /** Fills a range with random values */ + template + void generate(Iter first, Iter last) + { + for(; first != last; ++first) { + *first = (*this)(); + } + } -#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE +#ifndef BOOST_NO_LONG_LONG + /** + * Advances the state of the generator by @c z steps. Equivalent to + * + * @code + * for(unsigned long long i = 0; i < z; ++i) { + * gen(); + * } + * @endcode + */ + void discard(boost::ulong_long_type z) + { + for(boost::ulong_long_type j = 0; j < z; ++j) { + (*this)(); + } + } +#endif #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS - template - friend std::basic_ostream& - operator<<(std::basic_ostream& os, const mersenne_twister& mt) - { - for(int j = 0; j < mt.state_size; ++j) - os << mt.compute(j) << " "; - return os; - } - - template - friend std::basic_istream& - operator>>(std::basic_istream& is, mersenne_twister& mt) - { - for(int j = 0; j < mt.state_size; ++j) - is >> mt.x[j] >> std::ws; - // MSVC (up to 7.1) and Borland (up to 5.64) don't handle the template - // value parameter "n" available from the class template scope, so use - // the static constant with the same value - mt.i = mt.state_size; - return is; - } + /** Writes a mersenne_twister_engine to a @c std::ostream */ + template + friend std::basic_ostream& + operator<<(std::basic_ostream& os, + const mersenne_twister_engine& mt) + { + for(std::size_t j = 0; j < mt.state_size; ++j) + os << mt.compute(j) << " "; + return os; + } + + /** Reads a mersenne_twister_engine from a @c std::istream */ + template + friend std::basic_istream& + operator>>(std::basic_istream& is, mersenne_twister_engine& mt) + { + for(std::size_t j = 0; j < mt.state_size; ++j) + is >> mt.x[j] >> std::ws; + // MSVC (up to 7.1) and Borland (up to 5.64) don't handle the template + // value parameter "n" available from the class template scope, so use + // the static constant with the same value + mt.i = mt.state_size; + return is; + } #endif - friend bool operator==(const mersenne_twister& x, const mersenne_twister& y) - { - for(int j = 0; j < state_size; ++j) - if(x.compute(j) != y.compute(j)) - return false; - return true; - } - - friend bool operator!=(const mersenne_twister& x, const mersenne_twister& y) - { return !(x == y); } -#else - // Use a member function; Streamable concept not supported. - bool operator==(const mersenne_twister& rhs) const - { - for(int j = 0; j < state_size; ++j) - if(compute(j) != rhs.compute(j)) - return false; - return true; - } - - bool operator!=(const mersenne_twister& rhs) const - { return !(*this == rhs); } -#endif + /** + * Returns true if the two generators are in the same state, + * and will thus produce identical sequences. + */ + friend bool operator==(const mersenne_twister_engine& x, const mersenne_twister_engine& y) + { + for(std::size_t j = 0; j < state_size; ++j) + if(x.compute(j) != y.compute(j)) + return false; + return true; + } + + /** + * Returns true if the two generators are in different states. + */ + friend bool operator!=(const mersenne_twister_engine& x, const mersenne_twister_engine& y) + { return !(x == y); } private: - /// \cond hide_private_members - // returns x(i-n+index), where index is in 0..n-1 - UIntType compute(unsigned int index) const - { - // equivalent to (i-n+index) % 2n, but doesn't produce negative numbers - return x[ (i + n + index) % (2*n) ]; - } - void twist(int block); - /// \endcond + /// \cond hide_private_members + // returns x(i-n+index), where index is in 0..n-1 + UIntType compute(unsigned int index) const + { + // equivalent to (i-n+index) % 2n, but doesn't produce negative numbers + return x[ (i + n + index) % (2*n) ]; + } + void twist(int block); + /// \endcond - // state representation: next output is o(x(i)) - // x[0] ... x[k] x[k+1] ... x[n-1] x[n] ... x[2*n-1] represents - // x(i-k) ... x(i) x(i+1) ... x(i-k+n-1) x(i-k-n) ... x[i(i-k-1)] - // The goal is to always have x(i-n) ... x(i-1) available for - // operator== and save/restore. + // state representation: next output is o(x(i)) + // x[0] ... x[k] x[k+1] ... x[n-1] x[n] ... x[2*n-1] represents + // x(i-k) ... x(i) x(i+1) ... x(i-k+n-1) x(i-k-n) ... x[i(i-k-1)] + // The goal is to always have x(i-n) ... x(i-1) available for + // operator== and save/restore. - UIntType x[2*n]; - int i; + UIntType x[2*n]; + std::size_t i; }; #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION // A definition is required even for integral static constants -template -const bool mersenne_twister::has_fixed_range; -template -const int mersenne_twister::state_size; -template -const int mersenne_twister::shift_size; -template -const int mersenne_twister::mask_bits; -template -const UIntType mersenne_twister::parameter_a; -template -const int mersenne_twister::output_u; -template -const int mersenne_twister::output_s; -template -const UIntType mersenne_twister::output_b; -template -const int mersenne_twister::output_t; -template -const UIntType mersenne_twister::output_c; -template -const int mersenne_twister::output_l; +#define BOOST_RANDOM_MT_DEFINE_CONSTANT(type, name) \ +template \ +const type mersenne_twister_engine::name +BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, word_size); +BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, state_size); +BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, shift_size); +BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, mask_bits); +BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, xor_mask); +BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, tempering_u); +BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, tempering_d); +BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, tempering_s); +BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, tempering_b); +BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, tempering_t); +BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, tempering_c); +BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, tempering_l); +BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, initialization_multiplier); +BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, default_seed); +BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, parameter_a); +BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, output_u ); +BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, output_s); +BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, output_b); +BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, output_t); +BOOST_RANDOM_MT_DEFINE_CONSTANT(UIntType, output_c); +BOOST_RANDOM_MT_DEFINE_CONSTANT(std::size_t, output_l); +BOOST_RANDOM_MT_DEFINE_CONSTANT(bool, has_fixed_range); +#undef BOOST_RANDOM_MT_DEFINE_CONSTANT #endif /// \cond hide_private_members -template -void mersenne_twister::twist(int block) +template +void +mersenne_twister_engine::twist(int block) { - const UIntType upper_mask = (~0u) << r; - const UIntType lower_mask = ~upper_mask; + const UIntType upper_mask = (~static_cast(0)) << r; + const UIntType lower_mask = ~upper_mask; - if(block == 0) { - for(int j = n; j < 2*n; j++) { - UIntType y = (x[j-n] & upper_mask) | (x[j-(n-1)] & lower_mask); - x[j] = x[j-(n-m)] ^ (y >> 1) ^ (y&1 ? a : 0); - } - } else if (block == 1) { - // split loop to avoid costly modulo operations - { // extra scope for MSVC brokenness w.r.t. for scope - for(int j = 0; j < n-m; j++) { - UIntType y = (x[j+n] & upper_mask) | (x[j+n+1] & lower_mask); - x[j] = x[j+n+m] ^ (y >> 1) ^ (y&1 ? a : 0); - } - } + if(block == 0) { + for(std::size_t j = n; j < 2*n; j++) { + UIntType y = (x[j-n] & upper_mask) | (x[j-(n-1)] & lower_mask); + x[j] = x[j-(n-m)] ^ (y >> 1) ^ (y&1 ? a : 0); + } + } else if (block == 1) { + // split loop to avoid costly modulo operations + { // extra scope for MSVC brokenness w.r.t. for scope + for(std::size_t j = 0; j < n-m; j++) { + UIntType y = (x[j+n] & upper_mask) | (x[j+n+1] & lower_mask); + x[j] = x[j+n+m] ^ (y >> 1) ^ (y&1 ? a : 0); + } + } - for(int j = n-m; j < n-1; j++) { - UIntType y = (x[j+n] & upper_mask) | (x[j+n+1] & lower_mask); - x[j] = x[j-(n-m)] ^ (y >> 1) ^ (y&1 ? a : 0); + for(std::size_t j = n-m; j < n-1; j++) { + UIntType y = (x[j+n] & upper_mask) | (x[j+n+1] & lower_mask); + x[j] = x[j-(n-m)] ^ (y >> 1) ^ (y&1 ? a : 0); + } + // last iteration + UIntType y = (x[2*n-1] & upper_mask) | (x[0] & lower_mask); + x[n-1] = x[m-1] ^ (y >> 1) ^ (y&1 ? a : 0); + i = 0; } - // last iteration - UIntType y = (x[2*n-1] & upper_mask) | (x[0] & lower_mask); - x[n-1] = x[m-1] ^ (y >> 1) ^ (y&1 ? a : 0); - i = 0; - } } /// \endcond -template -inline typename mersenne_twister::result_type -mersenne_twister::operator()() +template +inline typename mersenne_twister_engine::result_type +mersenne_twister_engine::operator()() { - if(i == n) - twist(0); - else if(i >= 2*n) - twist(1); - // Step 4 - UIntType z = x[i]; - ++i; - z ^= (z >> u); - z ^= ((z << s) & b); - z ^= ((z << t) & c); - z ^= (z >> l); - return z; + if(i == n) + twist(0); + else if(i >= 2*n) + twist(1); + // Step 4 + UIntType z = x[i]; + ++i; + z ^= ((z >> u) & d); + z ^= ((z << s) & b); + z ^= ((z << t) & c); + z ^= (z >> l); + return z; } +/** + * The specializations \mt11213b and \mt19937 are from + * + * @blockquote + * "Mersenne Twister: A 623-dimensionally equidistributed + * uniform pseudo-random number generator", Makoto Matsumoto + * and Takuji Nishimura, ACM Transactions on Modeling and + * Computer Simulation: Special Issue on Uniform Random Number + * Generation, Vol. 8, No. 1, January 1998, pp. 3-30. + * @endblockquote + */ +typedef random::mersenne_twister_engine mt11213b; + +/** + * The specializations \mt11213b and \mt19937 are from + * + * @blockquote + * "Mersenne Twister: A 623-dimensionally equidistributed + * uniform pseudo-random number generator", Makoto Matsumoto + * and Takuji Nishimura, ACM Transactions on Modeling and + * Computer Simulation: Special Issue on Uniform Random Number + * Generation, Vol. 8, No. 1, January 1998, pp. 3-30. + * @endblockquote + */ +typedef random::mersenne_twister_engine mt19937; + +#if !defined(BOOST_NO_INT64_T) && !defined(BOOST_NO_INTEGRAL_INT64_T) +typedef random::mersenne_twister_engine mt19937_64; +#endif + } // namespace random -/** - * The specializations \mt11213b and \mt19937 are from - * - * @blockquote - * "Mersenne Twister: A 623-dimensionally equidistributed - * uniform pseudo-random number generator", Makoto Matsumoto - * and Takuji Nishimura, ACM Transactions on Modeling and - * Computer Simulation: Special Issue on Uniform Random Number - * Generation, Vol. 8, No. 1, January 1998, pp. 3-30. - * @endblockquote - */ -typedef random::mersenne_twister mt11213b; - -/** - * The specializations \mt11213b and \mt19937 are from - * - * @blockquote - * "Mersenne Twister: A 623-dimensionally equidistributed - * uniform pseudo-random number generator", Makoto Matsumoto - * and Takuji Nishimura, ACM Transactions on Modeling and - * Computer Simulation: Special Issue on Uniform Random Number - * Generation, Vol. 8, No. 1, January 1998, pp. 3-30. - * @endblockquote - */ -typedef random::mersenne_twister mt19937; +using random::mt11213b; +using random::mt19937; +using random::mt19937_64; } // namespace boost +BOOST_RANDOM_PTR_HELPER_SPEC(boost::mt11213b) BOOST_RANDOM_PTR_HELPER_SPEC(boost::mt19937) +BOOST_RANDOM_PTR_HELPER_SPEC(boost::mt19937_64) #endif // BOOST_RANDOM_MERSENNE_TWISTER_HPP diff --git a/include/boost/random/seed_seq.hpp b/include/boost/random/seed_seq.hpp new file mode 100644 index 0000000..cadfb93 --- /dev/null +++ b/include/boost/random/seed_seq.hpp @@ -0,0 +1,100 @@ +/* boost random/seed_seq.hpp header file + * + * Copyright Steven Watanabe 2010 + * Distributed under the Boost Software License, Version 1.0. (See + * accompanying file LICENSE_1_0.txt or copy at + * http://www.boost.org/LICENSE_1_0.txt) + * + * See http://www.boost.org for most recent version including documentation. + * + * $Id$ + * + */ + +#ifndef BOOST_RANDOM_SEED_SEQ_HPP +#define BOOST_RANDOM_SEED_SEQ_HPP + +#include +#include +#include +#include +#include +#include +#include +#include + +#ifndef BOOST_NO_INITIALIZER_LISTS +#include +#endif + +namespace boost { +namespace random { + +class seed_seq { +public: + typedef boost::uint_least32_t result_type; + + seed_seq() {} +#ifndef BOOST_NO_INITIALIZER_LISTS + template + seed_seq(const std::initializer_list& il) : v(il.begin(), il.end()) {} +#endif + template + seed_seq(Iter first, Iter last) : v(first, last) {} + template + explicit seed_seq(const Range& range) + : v(boost::begin(range), boost::end(range)) {} + + template + void generate(Iter first, Iter last) const + { + typedef typename std::iterator_traits::value_type value_type; + std::fill(first, last, static_cast(0x8b8b8b8bu)); + std::size_t s = v.size(); + std::size_t n = last - first; + std::size_t t = + (n >= 623) ? 11 : + (n >= 68) ? 7 : + (n >= 39) ? 5 : + (n >= 7) ? 3 : + (n - 1)/2; + std::size_t p = (n - t) / 2; + std::size_t q = p + t; + std::size_t m = (std::max)(s+1, n); + value_type mask = 0xffffffffu; + for(std::size_t k = 0; k < m; ++k) { + value_type r1 = + *(first + k%n) ^ *(first + (k+p)%n) ^ *(first + (k+n-1)%n); + r1 = r1 ^ (r1 >> 27); + r1 = (r1 * 1664525u) & mask; + value_type r2 = r1 + + ((k == 0) ? s : + (k <= s) ? k % n + v[k - 1] : + (k % n)); + *(first + (k+p)%n) = (*(first + (k+p)%n) + r1) & mask; + *(first + (k+q)%n) = (*(first + (k+q)%n) + r2) & mask; + *(first + k%n) = r2; + } + for(std::size_t k = m; k < m + n; ++k) { + value_type r3 = + (*(first + k%n) + *(first + (k+p)%n) + *(first + (k+n-1)%n)) + & mask; + r3 = r3 ^ (r3 >> 27); + r3 = (r3 * 1566083941u) & mask; + value_type r4 = r3 - k%m; + *(first + (k+p)%n) ^= r4; + *(first + (k+q)%n) ^= r3; + *(first + k%n) = r4; + } + } + std::size_t size() const { return v.size(); } + template + void param(Iter out) { std::copy(v.begin(), v.end(), out); } +private: + std::vector v; +}; + +} +} + +#endif diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 4f61341..3d18617 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -25,6 +25,7 @@ local all-urngs = hellekalek1995 mt11213b mt19937 + mt19937_64 lagged_fibonacci lagged_fibonacci607 ranlux3 @@ -41,6 +42,8 @@ for urng in $(all-urngs) run instantiate.cpp : : : BOOST_RANDOM_URNG_TEST=$(urng) : test_$(urng) ; } +run test_seed_seq.cpp /boost//unit_test_framework ; + run test_binomial.cpp ; run test_binomial_distribution.cpp /boost//unit_test_framework ; run test_poisson.cpp ; diff --git a/test/instantiate.cpp b/test/instantiate.cpp index 1ba024f..2bf926e 100644 --- a/test/instantiate.cpp +++ b/test/instantiate.cpp @@ -8,10 +8,19 @@ * $Id$ */ +#include + #if defined(BOOST_MSVC) && BOOST_MSVC <= 1300 #pragma warning( disable : 4786 ) #endif +#ifdef BOOST_MSVC +// conversion warnings. These are unavoidable, since we're +// testing all possible implicit conversions. +#pragma warning( disable : 4244 ) +#pragma warning( disable : 4305 ) +#endif + #include #include #include diff --git a/test/test_seed_seq.cpp b/test/test_seed_seq.cpp new file mode 100644 index 0000000..79a0e1a --- /dev/null +++ b/test/test_seed_seq.cpp @@ -0,0 +1,113 @@ +/* boost test_seed_seq.cpp + * + * Copyright Steven Watanabe 2010 + * Distributed under the Boost Software License, Version 1.0. (See + * accompanying file LICENSE_1_0.txt or copy at + * http://www.boost.org/LICENSE_1_0.txt) + * + * $Id$ + */ + +#include +#include +#include +#include + +#define BOOST_TEST_MAIN +#include + +using boost::assign::list_of; + +BOOST_AUTO_TEST_CASE(test_seed_seq) { + boost::uint32_t expected_param[4] = { 2, 3, 4, 0xdeadbeaf }; + boost::uint32_t param[4] = { 2, 3, 4, 0xdeadbeaf }; + boost::uint32_t store32[10]; + boost::uint64_t store64[10]; + boost::uint32_t expected[10] = { + 2571534471u, + 852093106u, + 3994144234u, + 904246505u, + 3404579686u, + 1508446734u, + 3679033230u, + 2550178480u, + 4205639301u, + 268720383u + }; + std::fill_n(&store32[0], 10, 0); + std::fill_n(&store64[0], 10, 0); + boost::random::seed_seq seq; + seq.generate(&store32[0], &store32[0] + 10); + BOOST_CHECK_EQUAL_COLLECTIONS( + &store32[0], &store32[0] + 10, &expected[0], &expected[0] + 10); + seq.generate(&store64[0], &store64[0] + 10); + BOOST_CHECK_EQUAL_COLLECTIONS( + &store64[0], &store64[0] + 10, &expected[0], &expected[0] + 10); + BOOST_CHECK_EQUAL(seq.size(), 0u); + seq.param(¶m[0]); + BOOST_CHECK_EQUAL_COLLECTIONS( + ¶m[0], ¶m[0] + 4, &expected_param[0], &expected_param[0] + 4); + + boost::uint32_t expected_r[10] = { + 249299179u, + 1602610063u, + 745496312u, + 3444440002u, + 1462959439u, + 4120764957u, + 2551678045u, + 2909251680u, + 3652199668u, + 1330108271u + }; + + std::vector data = list_of(2)(3)(4); + + std::fill_n(&store32[0], 10, 0); + std::fill_n(&store64[0], 10, 0); + std::fill_n(¶m[0], 3, 0); + boost::random::seed_seq seq_r(data); + seq_r.generate(&store32[0], &store32[0] + 10); + BOOST_CHECK_EQUAL_COLLECTIONS( + &store32[0], &store32[0] + 10, &expected_r[0], &expected_r[0] + 10); + seq_r.generate(&store64[0], &store64[0] + 10); + BOOST_CHECK_EQUAL_COLLECTIONS( + &store64[0], &store64[0] + 10, &expected_r[0], &expected_r[0] + 10); + BOOST_CHECK_EQUAL(seq_r.size(), 3u); + seq_r.param(¶m[0]); + BOOST_CHECK_EQUAL_COLLECTIONS( + ¶m[0], ¶m[0] + 4, &expected_param[0], &expected_param[0] + 4); + + std::fill_n(&store32[0], 10, 0); + std::fill_n(&store64[0], 10, 0); + std::fill_n(¶m[0], 3, 0); + boost::random::seed_seq seq_it(data.begin(), data.end()); + seq_it.generate(&store32[0], &store32[0] + 10); + BOOST_CHECK_EQUAL_COLLECTIONS( + &store32[0], &store32[0] + 10, &expected_r[0], &expected_r[0] + 10); + seq_it.generate(&store64[0], &store64[0] + 10); + BOOST_CHECK_EQUAL_COLLECTIONS( + &store64[0], &store64[0] + 10, &expected_r[0], &expected_r[0] + 10); + BOOST_CHECK_EQUAL(seq_it.size(), 3u); + seq_it.param(¶m[0]); + BOOST_CHECK_EQUAL_COLLECTIONS( + ¶m[0], ¶m[0] + 4, &expected_param[0], &expected_param[0] + 4); + +#ifndef BOOST_NO_INITIALIZER_LISTS + std::fill_n(&store32[0], 10, 0); + std::fill_n(&store64[0], 10, 0); + std::fill_n(¶m[0], 3, 0); + boost::random::seed_seq seq_il = {2, 3, 4}; + seq_il.generate(&store32[0], &store32[0] + 10); + BOOST_CHECK_EQUAL_COLLECTIONS( + &store32[0], &store32[0] + 10, &expected_r[0], &expected_r[0] + 10); + seq_il.generate(&store64[0], &store64[0] + 10); + BOOST_CHECK_EQUAL_COLLECTIONS( + &store64[0], &store64[0] + 10, &expected_r[0], &expected_r[0] + 10); + BOOST_CHECK_EQUAL(seq_il.size(), 3u); + seq_il.param(¶m[0]); + BOOST_CHECK_EQUAL_COLLECTIONS( + ¶m[0], ¶m[0] + 4, &expected_param[0], &expected_param[0] + 4); +#endif +} diff --git a/test/validate.cpp b/test/validate.cpp index 2a38cf2..a4bd728 100644 --- a/test/validate.cpp +++ b/test/validate.cpp @@ -49,6 +49,10 @@ bool check_(unsigned long x, const boost::mt11213b&) { return x == 3809585648U; // validation by experiment from mt19937.c bool check_(unsigned long x, const boost::mt19937&) { return x == 4123659995U; } +#if !defined(BOOST_NO_INT64_T) && !defined(BOOST_NO_INTEGRAL_INT64_T) +bool check_(boost::uint64_t x, const boost::mt19937_64&) { return x == UINT64_C(9981545732273789042); } +#endif + // validation values from the publications bool check_(int x, const boost::minstd_rand0&) { return x == 1043618065; } @@ -113,6 +117,9 @@ void validate_all() validate("ecuyer combined", ecuyer1988()); validate("mt11213b", mt11213b()); validate("mt19937", mt19937()); +#if !defined(BOOST_NO_INT64_T) && !defined(BOOST_NO_INTEGRAL_INT64_T) + validate("mt19937_64", mt19937_64()); +#endif validate("kreutzer1986", kreutzer1986()); validate("ranlux3", ranlux3()); validate("ranlux4", ranlux4());