Swiftdiff options
Diffstat (limited to '3rdParty/Boost/src/boost/random')
13 files changed, 1875 insertions, 976 deletions
diff --git a/3rdParty/Boost/src/boost/random/detail/const_mod.hpp b/3rdParty/Boost/src/boost/random/detail/const_mod.hpp index e0a8839..9778f55 100644 --- a/3rdParty/Boost/src/boost/random/detail/const_mod.hpp +++ b/3rdParty/Boost/src/boost/random/detail/const_mod.hpp @@ -7,7 +7,7 @@ * * See http://www.boost.org for most recent version including documentation. * - * $Id: const_mod.hpp 58649 2010-01-02 21:23:17Z steven_watanabe $ + * $Id: const_mod.hpp 71018 2011-04-05 21:27:52Z steven_watanabe $ * * Revision history * 2001-02-18 moved to individual header files @@ -16,113 +16,101 @@ #ifndef BOOST_RANDOM_CONST_MOD_HPP #define BOOST_RANDOM_CONST_MOD_HPP -#include <cassert> +#include <boost/assert.hpp> #include <boost/static_assert.hpp> -#include <boost/cstdint.hpp> #include <boost/integer_traits.hpp> -#include <boost/detail/workaround.hpp> +#include <boost/type_traits/make_unsigned.hpp> +#include <boost/random/detail/large_arithmetic.hpp> #include <boost/random/detail/disable_warnings.hpp> namespace boost { namespace random { -/* - * Some random number generators require modular arithmetic. Put - * everything we need here. - * IntType must be an integral type. - */ - -namespace detail { - - template<bool is_signed> - struct do_add - { }; - - template<> - struct do_add<true> - { - template<class IntType> - static IntType add(IntType m, IntType x, IntType c) - { - if (x < m - c) - return x + c; - else - return x - (m-c); - } - }; - - template<> - struct do_add<false> - { - template<class IntType> - static IntType add(IntType, IntType, IntType) - { - // difficult - assert(!"const_mod::add with c too large"); - return 0; - } - }; -} // namespace detail - -#if !(defined(__BORLANDC__) && (__BORLANDC__ == 0x560)) - template<class IntType, IntType m> class const_mod { public: + static IntType apply(IntType x) + { + if(((unsigned_m() - 1) & unsigned_m()) == 0) + return (unsigned_type(x)) & (unsigned_m() - 1); + else { + IntType supress_warnings = (m == 0); + BOOST_ASSERT(supress_warnings == 0); + return x % (m + supress_warnings); + } + } + static IntType add(IntType x, IntType c) { - if(c == 0) + if(((unsigned_m() - 1) & unsigned_m()) == 0) + return (unsigned_type(x) + unsigned_type(c)) & (unsigned_m() - 1); + else if(c == 0) return x; - else if(c <= traits::const_max - m) // i.e. m+c < max - return add_small(x, c); + else if(x < m - c) + return x + c; else - return detail::do_add<traits::is_signed>::add(m, x, c); + return x - (m - c); } static IntType mult(IntType a, IntType x) { - if(a == 1) + if(((unsigned_m() - 1) & unsigned_m()) == 0) + return unsigned_type(a) * unsigned_type(x) & (unsigned_m() - 1); + else if(a == 0) + return 0; + else if(a == 1) return x; else if(m <= traits::const_max/a) // i.e. a*m <= max return mult_small(a, x); else if(traits::is_signed && (m%a < m/a)) return mult_schrage(a, x); - else { - // difficult - assert(!"const_mod::mult with a too large"); - return 0; - } + else + return mult_general(a, x); } static IntType mult_add(IntType a, IntType x, IntType c) { - if(m <= (traits::const_max-c)/a) // i.e. a*m+c <= max - return (a*x+c) % m; - else + if(((unsigned_m() - 1) & unsigned_m()) == 0) + return (unsigned_type(a) * unsigned_type(x) + unsigned_type(c)) & (unsigned_m() - 1); + else if(a == 0) + return c; + else if(m <= (traits::const_max-c)/a) { // i.e. a*m+c <= max + IntType supress_warnings = (m == 0); + BOOST_ASSERT(supress_warnings == 0); + return (a*x+c) % (m + supress_warnings); + } else return add(mult(a, x), c); } + static IntType pow(IntType a, boost::uintmax_t exponent) + { + IntType result = 1; + while(exponent != 0) { + if(exponent % 2 == 1) { + result = mult(result, a); + } + a = mult(a, a); + exponent /= 2; + } + return result; + } + static IntType invert(IntType x) - { return x == 0 ? 0 : invert_euclidian(x); } + { return x == 0 ? 0 : (m == 0? invert_euclidian0(x) : invert_euclidian(x)); } private: typedef integer_traits<IntType> traits; + typedef typename make_unsigned<IntType>::type unsigned_type; const_mod(); // don't instantiate - static IntType add_small(IntType x, IntType c) - { - x += c; - if(x >= m) - x -= m; - return x; - } - static IntType mult_small(IntType a, IntType x) { - return a*x % m; + IntType supress_warnings = (m == 0); + BOOST_ASSERT(supress_warnings == 0); + return a*x % (m + supress_warnings); } static IntType mult_schrage(IntType a, IntType value) @@ -130,231 +118,96 @@ private: const IntType q = m / a; const IntType r = m % a; - assert(r < q); // check that overflow cannot happen + BOOST_ASSERT(r < q); // check that overflow cannot happen - value = a*(value%q) - r*(value/q); - // An optimizer bug in the SGI MIPSpro 7.3.1.x compiler requires this - // convoluted formulation of the loop (Synge Todo) - for(;;) { - if (value > 0) - break; - value += m; + return sub(a*(value%q), r*(value/q)); + } + + static IntType mult_general(IntType a, IntType b) + { + IntType suppress_warnings = (m == 0); + BOOST_ASSERT(suppress_warnings == 0); + IntType modulus = m + suppress_warnings; + BOOST_ASSERT(modulus == m); + if(::boost::uintmax_t(modulus) <= + (::std::numeric_limits< ::boost::uintmax_t>::max)() / modulus) + { + return static_cast<IntType>(boost::uintmax_t(a) * b % modulus); + } else { + return static_cast<IntType>(detail::mulmod(a, b, modulus)); } - return value; + } + + static IntType sub(IntType a, IntType b) + { + if(a < b) + return m - (b - a); + else + return a - b; + } + + static unsigned_type unsigned_m() + { + if(m == 0) { + return unsigned_type((std::numeric_limits<IntType>::max)()) + 1; + } else { + return unsigned_type(m); + } } // invert c in the finite field (mod m) (m must be prime) static IntType invert_euclidian(IntType c) { // we are interested in the gcd factor for c, because this is our inverse - BOOST_STATIC_ASSERT(m > 0); -#if BOOST_WORKAROUND(__MWERKS__, BOOST_TESTED_AT(0x3003)) - assert(boost::integer_traits<IntType>::is_signed); -#elif !defined(BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS) - BOOST_STATIC_ASSERT(boost::integer_traits<IntType>::is_signed); -#endif - assert(c > 0); + BOOST_ASSERT(c > 0); IntType l1 = 0; IntType l2 = 1; IntType n = c; IntType p = m; for(;;) { IntType q = p / n; - l1 -= q * l2; // this requires a signed IntType! + l1 += q * l2; p -= q * n; if(p == 0) - return (l2 < 1 ? l2 + m : l2); + return l2; IntType q2 = n / p; - l2 -= q2 * l1; + l2 += q2 * l1; n -= q2 * p; if(n == 0) - return (l1 < 1 ? l1 + m : l1); + return m - l1; } } -}; - -// The modulus is exactly the word size: rely on machine overflow handling. -// Due to a GCC bug, we cannot partially specialize in the presence of -// template value parameters. -template<> -class const_mod<unsigned int, 0> -{ - typedef unsigned int IntType; -public: - static IntType add(IntType x, IntType c) { return x+c; } - static IntType mult(IntType a, IntType x) { return a*x; } - static IntType mult_add(IntType a, IntType x, IntType c) { return a*x+c; } - - // m is not prime, thus invert is not useful -private: // don't instantiate - const_mod(); -}; - -template<> -class const_mod<unsigned long, 0> -{ - typedef unsigned long IntType; -public: - static IntType add(IntType x, IntType c) { return x+c; } - static IntType mult(IntType a, IntType x) { return a*x; } - static IntType mult_add(IntType a, IntType x, IntType c) { return a*x+c; } - - // m is not prime, thus invert is not useful -private: // don't instantiate - const_mod(); -}; - -// the modulus is some power of 2: rely partly on machine overflow handling -// we only specialize for rand48 at the moment -#ifndef BOOST_NO_INT64_T -template<> -class const_mod<uint64_t, uint64_t(1) << 48> -{ - typedef uint64_t IntType; -public: - static IntType add(IntType x, IntType c) { return c == 0 ? x : mod(x+c); } - static IntType mult(IntType a, IntType x) { return mod(a*x); } - static IntType mult_add(IntType a, IntType x, IntType c) - { return mod(a*x+c); } - static IntType mod(IntType x) { return x &= ((uint64_t(1) << 48)-1); } - - // m is not prime, thus invert is not useful -private: // don't instantiate - const_mod(); -}; -#endif /* !BOOST_NO_INT64_T */ - -#else - -// -// for some reason Borland C++ Builder 6 has problems with -// the full specialisations of const_mod, define a generic version -// instead, the compiler will optimise away the const-if statements: -// -template<class IntType, IntType m> -class const_mod -{ -public: - static IntType add(IntType x, IntType c) - { - if(0 == m) - { - return x+c; - } - else - { - if(c == 0) - return x; - else if(c <= traits::const_max - m) // i.e. m+c < max - return add_small(x, c); - else - return detail::do_add<traits::is_signed>::add(m, x, c); - } - } - - static IntType mult(IntType a, IntType x) - { - if(x == 0) - { - return a*x; - } - else - { - if(a == 1) - return x; - else if(m <= traits::const_max/a) // i.e. a*m <= max - return mult_small(a, x); - else if(traits::is_signed && (m%a < m/a)) - return mult_schrage(a, x); - else { - // difficult - assert(!"const_mod::mult with a too large"); - return 0; - } - } - } - - static IntType mult_add(IntType a, IntType x, IntType c) - { - if(m == 0) - { - return a*x+c; - } - else - { - if(m <= (traits::const_max-c)/a) // i.e. a*m+c <= max - return (a*x+c) % m; - else - return add(mult(a, x), c); - } - } - - static IntType invert(IntType x) - { return x == 0 ? 0 : invert_euclidian(x); } - -private: - typedef integer_traits<IntType> traits; - - const_mod(); // don't instantiate - - static IntType add_small(IntType x, IntType c) - { - x += c; - if(x >= m) - x -= m; - return x; - } - - static IntType mult_small(IntType a, IntType x) - { - return a*x % m; - } - - static IntType mult_schrage(IntType a, IntType value) - { - const IntType q = m / a; - const IntType r = m % a; - - assert(r < q); // check that overflow cannot happen - - value = a*(value%q) - r*(value/q); - while(value <= 0) - value += m; - return value; - } - - // invert c in the finite field (mod m) (m must be prime) - static IntType invert_euclidian(IntType c) + // invert c in the finite field (mod m) (c must be relatively prime to m) + static IntType invert_euclidian0(IntType c) { // we are interested in the gcd factor for c, because this is our inverse - BOOST_STATIC_ASSERT(m > 0); -#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS - BOOST_STATIC_ASSERT(boost::integer_traits<IntType>::is_signed); -#endif - assert(c > 0); + BOOST_ASSERT(c > 0); + if(c == 1) return 1; IntType l1 = 0; IntType l2 = 1; IntType n = c; IntType p = m; + IntType max = (std::numeric_limits<IntType>::max)(); + IntType q = max / n; + BOOST_ASSERT(max % n != n - 1 && "c must be relatively prime to m."); + l1 += q * l2; + p = max - q * n + 1; for(;;) { - IntType q = p / n; - l1 -= q * l2; // this requires a signed IntType! - p -= q * n; if(p == 0) - return (l2 < 1 ? l2 + m : l2); + return l2; IntType q2 = n / p; - l2 -= q2 * l1; + l2 += q2 * l1; n -= q2 * p; if(n == 0) - return (l1 < 1 ? l1 + m : l1); + return m - l1; + q = p / n; + l1 += q * l2; + p -= q * n; } } }; - -#endif - } // namespace random } // namespace boost diff --git a/3rdParty/Boost/src/boost/random/detail/generator_bits.hpp b/3rdParty/Boost/src/boost/random/detail/generator_bits.hpp new file mode 100644 index 0000000..44b4248 --- /dev/null +++ b/3rdParty/Boost/src/boost/random/detail/generator_bits.hpp @@ -0,0 +1,36 @@ +/* boost random/detail/generator_bits.hpp header file + * + * Copyright Steven Watanabe 2011 + * 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: generator_bits.hpp 72951 2011-07-07 04:57:37Z steven_watanabe $ + * + */ + +#ifndef BOOST_RANDOM_DETAIL_GENERATOR_BITS_HPP +#define BOOST_RANDOM_DETAIL_GENERATOR_BITS_HPP + +#include <boost/limits.hpp> + +namespace boost { +namespace random { +namespace detail { + +// This is a temporary measure that retains backwards +// compatibility. +template<class URNG> +struct generator_bits { + static std::size_t value() { + return std::numeric_limits<typename URNG::result_type>::digits; + } +}; + +} // namespace detail +} // namespace random +} // namespace boost + +#endif // BOOST_RANDOM_DETAIL_GENERATOR_BITS_HPP diff --git a/3rdParty/Boost/src/boost/random/detail/generator_seed_seq.hpp b/3rdParty/Boost/src/boost/random/detail/generator_seed_seq.hpp new file mode 100644 index 0000000..6aaf98f --- /dev/null +++ b/3rdParty/Boost/src/boost/random/detail/generator_seed_seq.hpp @@ -0,0 +1,40 @@ +/* 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) + * + * See http://www.boost.org for most recent version including documentation. + * + * $Id: generator_seed_seq.hpp 71018 2011-04-05 21:27:52Z steven_watanabe $ + * + */ + +#ifndef BOOST_RANDOM_DETAIL_GENERATOR_SEED_SEQ_HPP_INCLUDED +#define BOOST_RANDOM_DETAIL_GENERATOR_SEED_SEQ_HPP_INCLUDED + +namespace boost { +namespace random { +namespace detail { + +template<class Generator> +class generator_seed_seq { +public: + generator_seed_seq(Generator& g) : gen(&g) {} + template<class It> + void generate(It first, It last) { + for(; first != last; ++first) { + *first = (*gen)(); + } + } +private: + Generator* gen; +}; + +} +} +} + +#endif diff --git a/3rdParty/Boost/src/boost/random/detail/integer_log2.hpp b/3rdParty/Boost/src/boost/random/detail/integer_log2.hpp new file mode 100644 index 0000000..3770fc0 --- /dev/null +++ b/3rdParty/Boost/src/boost/random/detail/integer_log2.hpp @@ -0,0 +1,72 @@ +/* boost random/detail/integer_log2.hpp header file + * + * Copyright Steven Watanabe 2011 + * 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: integer_log2.hpp 72861 2011-07-02 20:26:19Z danieljames $ + * + */ + +#ifndef BOOST_RANDOM_DETAIL_INTEGER_LOG2_HPP +#define BOOST_RANDOM_DETAIL_INTEGER_LOG2_HPP + +#include <boost/config.hpp> +#include <boost/limits.hpp> +#include <boost/pending/integer_log2.hpp> + +namespace boost { +namespace random { +namespace detail { + +// Daniel James: Disabled use of constexpr because integer_log2_impl is not a +// valid constexpr. +#if 0 && !defined(BOOST_NO_CONSTEXPR) +#define BOOST_RANDOM_DETAIL_CONSTEXPR constexpr +#elif defined(BOOST_MSVC) +#define BOOST_RANDOM_DETAIL_CONSTEXPR __forceinline +#elif defined(__GNUC__) && __GNUC__ >= 4 +#define BOOST_RANDOM_DETAIL_CONSTEXPR __attribute__((const)) __attribute__((always_inline)) +#else +#define BOOST_RANDOM_DETAIL_CONSTEXPR inline +#endif + +template<int Shift> +struct integer_log2_impl +{ + template<class T> + BOOST_RANDOM_DETAIL_CONSTEXPR static int apply(T t, int accum) + { + int update = ((t >> Shift) != 0) * Shift; + return integer_log2_impl<Shift / 2>::apply(t >> update, accum + update); + } +}; + +template<> +struct integer_log2_impl<1> +{ + template<class T> + BOOST_RANDOM_DETAIL_CONSTEXPR static int apply(T t, int accum) + { + return int(t >> 1) + accum; + } +}; + +template<class T> +BOOST_RANDOM_DETAIL_CONSTEXPR int integer_log2(T t) +{ + return integer_log2_impl< + ::boost::detail::max_pow2_less< + ::std::numeric_limits<T>::digits, 4 + >::value + >::apply(t, 0); +} + +} // namespace detail +} // namespace random +} // namespace boost + +#endif // BOOST_RANDOM_DETAIL_INTEGER_LOG2_HPP diff --git a/3rdParty/Boost/src/boost/random/detail/large_arithmetic.hpp b/3rdParty/Boost/src/boost/random/detail/large_arithmetic.hpp new file mode 100644 index 0000000..24177dc --- /dev/null +++ b/3rdParty/Boost/src/boost/random/detail/large_arithmetic.hpp @@ -0,0 +1,122 @@ +/* boost random/detail/large_arithmetic.hpp header file + * + * Copyright Steven Watanabe 2011 + * 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: large_arithmetic.hpp 71018 2011-04-05 21:27:52Z steven_watanabe $ + */ + +#ifndef BOOST_RANDOM_DETAIL_LARGE_ARITHMETIC_HPP +#define BOOST_RANDOM_DETAIL_LARGE_ARITHMETIC_HPP + +#include <boost/cstdint.hpp> +#include <boost/integer.hpp> +#include <boost/limits.hpp> +#include <boost/random/detail/integer_log2.hpp> + +#include <boost/random/detail/disable_warnings.hpp> + +namespace boost { +namespace random { +namespace detail { + +struct div_t { + boost::uintmax_t quotient; + boost::uintmax_t remainder; +}; + +inline div_t muldivmod(boost::uintmax_t a, boost::uintmax_t b, boost::uintmax_t m) +{ + static const int bits = + ::std::numeric_limits< ::boost::uintmax_t>::digits / 2; + static const ::boost::uintmax_t mask = (::boost::uintmax_t(1) << bits) - 1; + typedef ::boost::uint_t<bits>::fast digit_t; + + int shift = std::numeric_limits< ::boost::uintmax_t>::digits - 1 + - detail::integer_log2(m); + + a <<= shift; + m <<= shift; + + digit_t product[4] = { 0, 0, 0, 0 }; + digit_t a_[2] = { digit_t(a & mask), digit_t((a >> bits) & mask) }; + digit_t b_[2] = { digit_t(b & mask), digit_t((b >> bits) & mask) }; + digit_t m_[2] = { digit_t(m & mask), digit_t((m >> bits) & mask) }; + + // multiply a * b + for(int i = 0; i < 2; ++i) { + digit_t carry = 0; + for(int j = 0; j < 2; ++j) { + ::boost::uint64_t temp = ::boost::uintmax_t(a_[i]) * b_[j] + + carry + product[i + j]; + product[i + j] = digit_t(temp & mask); + carry = digit_t(temp >> bits); + } + if(carry != 0) { + product[i + 2] += carry; + } + } + + digit_t quotient[2]; + + if(m == 0) { + div_t result = { + ((::boost::uintmax_t(product[3]) << bits) | product[2]), + ((::boost::uintmax_t(product[1]) << bits) | product[0]) >> shift, + }; + return result; + } + + // divide product / m + for(int i = 3; i >= 2; --i) { + ::boost::uintmax_t temp = + ::boost::uintmax_t(product[i]) << bits | product[i - 1]; + + digit_t q = digit_t((product[i] == m_[1]) ? mask : temp / m_[1]); + + ::boost::uintmax_t rem = + ((temp - ::boost::uintmax_t(q) * m_[1]) << bits) + product[i - 2]; + + ::boost::uintmax_t diff = m_[0] * ::boost::uintmax_t(q); + + int error = 0; + if(diff > rem) { + if(diff - rem > m) { + error = 2; + } else { + error = 1; + } + } + q -= error; + rem = rem + error * m - diff; + + quotient[i - 2] = q; + product[i] = 0; + product[i-1] = (rem >> bits) & mask; + product[i-2] = rem & mask; + } + + div_t result = { + ((::boost::uintmax_t(quotient[1]) << bits) | quotient[0]), + ((::boost::uintmax_t(product[1]) << bits) | product[0]) >> shift, + }; + return result; +} + +inline boost::uintmax_t muldiv(boost::uintmax_t a, boost::uintmax_t b, boost::uintmax_t m) +{ return detail::muldivmod(a, b, m).quotient; } + +inline boost::uintmax_t mulmod(boost::uintmax_t a, boost::uintmax_t b, boost::uintmax_t m) +{ return detail::muldivmod(a, b, m).remainder; } + +} // namespace detail +} // namespace random +} // namespace boost + +#include <boost/random/detail/enable_warnings.hpp> + +#endif // BOOST_RANDOM_DETAIL_LARGE_ARITHMETIC_HPP diff --git a/3rdParty/Boost/src/boost/random/detail/operators.hpp b/3rdParty/Boost/src/boost/random/detail/operators.hpp new file mode 100644 index 0000000..f27839a --- /dev/null +++ b/3rdParty/Boost/src/boost/random/detail/operators.hpp @@ -0,0 +1,84 @@ +/* boost random/detail/operators.hpp header file + * + * Copyright Steven Watanabe 2010-2011 + * 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: operators.hpp 71018 2011-04-05 21:27:52Z steven_watanabe $ + */ + +#ifndef BOOST_RANDOM_DETAIL_OPERATORS_HPP +#define BOOST_RANDOM_DETAIL_OPERATORS_HPP + +#include <boost/random/detail/config.hpp> +#include <boost/detail/workaround.hpp> + +#if BOOST_WORKAROUND(BOOST_MSVC, <= 1310) \ + || BOOST_WORKAROUND(__SUNPRO_CC, BOOST_TESTED_AT(0x5100)) + +#define BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, T, t) \ + template<class CharT, class Traits> \ + friend std::basic_ostream<CharT,Traits>& \ + operator<<(std::basic_ostream<CharT,Traits>& os, const T& t) { \ + t.print(os, t); \ + return os; \ + } \ + template<class CharT, class Traits> \ + static std::basic_ostream<CharT,Traits>& \ + print(std::basic_ostream<CharT,Traits>& os, const T& t) + +#define BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, T, t) \ + template<class CharT, class Traits> \ + friend std::basic_istream<CharT,Traits>& \ + operator>>(std::basic_istream<CharT,Traits>& is, T& t) { \ + t.read(is, t); \ + return is; \ + } \ + template<class CharT, class Traits> \ + static std::basic_istream<CharT,Traits>& \ + read(std::basic_istream<CharT,Traits>& is, T& t) + +#endif + +#if defined(__BORLANDC__) + +#define BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(T, lhs, rhs) \ + bool operator==(const T& rhs) const \ + { return T::is_equal(*this, rhs); } \ + static bool is_equal(const T& lhs, const T& rhs) + +#define BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(T) \ + bool operator!=(const T& rhs) const \ + { return !T::is_equal(*this, rhs); } + +#endif + +#ifndef BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR +#define BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, T, t) \ + template<class CharT, class Traits> \ + friend std::basic_ostream<CharT,Traits>& \ + operator<<(std::basic_ostream<CharT,Traits>& os, const T& t) +#endif + +#ifndef BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR +#define BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, T, t) \ + template<class CharT, class Traits> \ + friend std::basic_istream<CharT,Traits>& \ + operator>>(std::basic_istream<CharT,Traits>& is, T& t) +#endif + +#ifndef BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR +#define BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(T, lhs, rhs) \ + friend bool operator==(const T& lhs, const T& rhs) +#endif + +#ifndef BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR +#define BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(T) \ + friend bool operator!=(const T& lhs, const T& rhs) \ + { return !(lhs == rhs); } +#endif + +#endif diff --git a/3rdParty/Boost/src/boost/random/detail/seed.hpp b/3rdParty/Boost/src/boost/random/detail/seed.hpp index 48cc17e..979db29 100644 --- a/3rdParty/Boost/src/boost/random/detail/seed.hpp +++ b/3rdParty/Boost/src/boost/random/detail/seed.hpp @@ -7,7 +7,7 @@ * * See http://www.boost.org for most recent version including documentation. * - * $Id: seed.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * $Id: seed.hpp 71018 2011-04-05 21:27:52Z steven_watanabe $ */ #ifndef BOOST_RANDOM_DETAIL_SEED_HPP @@ -43,12 +43,19 @@ struct disable_constructor<Engine, Engine> {}; template<class Generator> \ void seed(Generator& gen, typename ::boost::random::detail::disable_seed<Generator>::type* = 0) +#define BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(Self, SeedSeq, seq) \ + template<class SeedSeq> \ + explicit Self(SeedSeq& seq, typename ::boost::random::detail::disable_constructor<Self, SeedSeq>::type* = 0) + +#define BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(Self, SeedSeq, seq) \ + template<class SeedSeq> \ + void seed(SeedSeq& seq, typename ::boost::random::detail::disable_seed<SeedSeq>::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<Engine, Engine> {}; template<class Generator>\ 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<class SeedSeq> \ + explicit Self(SeedSeq& seq) { \ + boost_random_constructor_impl(seq, ::boost::is_arithmetic<SeedSeq>());\ + } \ + template<class SeedSeq> \ + void boost_random_constructor_impl(SeedSeq& seq, ::boost::mpl::false_) + +#define BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(Self, SeedSeq, seq) \ + template<class SeedSeq> \ + void seed(SeedSeq& seq) { \ + boost_random_seed_impl(seq, ::boost::is_arithmetic<SeedSeq>()); \ + } \ + template<class SeedSeq> \ + 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/3rdParty/Boost/src/boost/random/detail/seed_impl.hpp b/3rdParty/Boost/src/boost/random/detail/seed_impl.hpp new file mode 100644 index 0000000..e044d45 --- /dev/null +++ b/3rdParty/Boost/src/boost/random/detail/seed_impl.hpp @@ -0,0 +1,397 @@ +/* boost random/detail/seed.hpp header file + * + * Copyright Steven Watanabe 2009 + * 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: seed_impl.hpp 72951 2011-07-07 04:57:37Z steven_watanabe $ + */ + +#ifndef BOOST_RANDOM_DETAIL_SEED_IMPL_HPP +#define BOOST_RANDOM_DETAIL_SEED_IMPL_HPP + +#include <stdexcept> +#include <boost/cstdint.hpp> +#include <boost/config/no_tr1/cmath.hpp> +#include <boost/integer/integer_mask.hpp> +#include <boost/integer/static_log2.hpp> +#include <boost/type_traits/is_signed.hpp> +#include <boost/type_traits/is_integral.hpp> +#include <boost/type_traits/make_unsigned.hpp> +#include <boost/mpl/bool.hpp> +#include <boost/mpl/if.hpp> +#include <boost/mpl/int.hpp> +#include <boost/random/detail/const_mod.hpp> +#include <boost/random/detail/integer_log2.hpp> +#include <boost/random/detail/signed_unsigned_tools.hpp> +#include <boost/random/detail/generator_bits.hpp> + +#include <boost/random/detail/disable_warnings.hpp> + +namespace boost { +namespace random { +namespace detail { + +// finds the seed type of an engine, given its +// result_type. If the result_type is integral +// the seed type is the same. If the result_type +// is floating point, the seed type is uint32_t +template<class T> +struct seed_type +{ + typedef typename boost::mpl::if_<boost::is_integral<T>, + T, + boost::uint32_t + >::type type; +}; + +template<int N> +struct const_pow_impl +{ + template<class T> + static T call(T arg, int n, T result) + { + return const_pow_impl<N / 2>::call(arg * arg, n / 2, + n%2 == 0? result : result * arg); + } +}; + +template<> +struct const_pow_impl<0> +{ + template<class T> + static T call(T, int, T result) + { + return result; + } +}; + +// requires N is an upper bound on n +template<int N, class T> +inline T const_pow(T arg, int n) { return const_pow_impl<N>::call(arg, n, T(1)); } + +template<class T> +inline T pow2(int n) +{ + typedef unsigned int_type; + const int max_bits = std::numeric_limits<int_type>::digits; + T multiplier = T(int_type(1) << (max_bits - 1)) * 2; + return (int_type(1) << (n % max_bits)) * + const_pow<std::numeric_limits<T>::digits / max_bits>(multiplier, n / max_bits); +} + +template<class Engine, class Iter> +void generate_from_real(Engine& eng, Iter begin, Iter end) +{ + using std::fmod; + typedef typename Engine::result_type RealType; + const int Bits = detail::generator_bits<Engine>::value(); + int remaining_bits = 0; + boost::uint_least32_t saved_bits = 0; + RealType multiplier = pow2<RealType>( Bits); + RealType mult32 = RealType(4294967296.0); // 2^32 + while(true) { + RealType val = eng() * multiplier; + int available_bits = Bits; + // Make sure the compiler can optimize this out + // if it isn't possible. + if(Bits < 32 && available_bits < 32 - remaining_bits) { + saved_bits |= boost::uint_least32_t(val) << remaining_bits; + remaining_bits += Bits; + } else { + // If Bits < 32, then remaining_bits != 0, since + // if remaining_bits == 0, available_bits < 32 - 0, + // and we won't get here to begin with. + if(Bits < 32 || remaining_bits != 0) { + boost::uint_least32_t divisor = + (boost::uint_least32_t(1) << (32 - remaining_bits)); + boost::uint_least32_t extra_bits = boost::uint_least32_t(fmod(val, mult32)) & (divisor - 1); + val = val / divisor; + *begin++ = saved_bits | (extra_bits << remaining_bits); + if(begin == end) return; + available_bits -= 32 - remaining_bits; + remaining_bits = 0; + } + // If Bits < 32 we should never enter this loop + if(Bits >= 32) { + for(; available_bits >= 32; available_bits -= 32) { + boost::uint_least32_t word = boost::uint_least32_t(fmod(val, mult32)); + val /= mult32; + *begin++ = word; + if(begin == end) return; + } + } + remaining_bits = available_bits; + saved_bits = static_cast<boost::uint_least32_t>(val); + } + } +} + +template<class Engine, class Iter> +void generate_from_int(Engine& eng, Iter begin, Iter end) +{ + typedef typename Engine::result_type IntType; + typedef typename boost::make_unsigned<IntType>::type unsigned_type; + int remaining_bits = 0; + boost::uint_least32_t saved_bits = 0; + unsigned_type range = boost::random::detail::subtract<IntType>()((eng.max)(), (eng.min)()); + + int bits = + (range == (std::numeric_limits<unsigned_type>::max)()) ? + std::numeric_limits<unsigned_type>::digits : + detail::integer_log2(range + 1); + + { + int discarded_bits = detail::integer_log2(bits); + unsigned_type excess = (range + 1) >> (bits - discarded_bits); + if(excess != 0) { + int extra_bits = detail::integer_log2((excess - 1) ^ excess); + bits = bits - discarded_bits + extra_bits; + } + } + + unsigned_type mask = (static_cast<unsigned_type>(2) << (bits - 1)) - 1; + unsigned_type limit = ((range + 1) & ~mask) - 1; + + while(true) { + unsigned_type val; + do { + val = boost::random::detail::subtract<IntType>()(eng(), (eng.min)()); + } while(limit != range && val > limit); + val &= mask; + int available_bits = bits; + if(available_bits == 32) { + *begin++ = static_cast<boost::uint_least32_t>(val) & 0xFFFFFFFFu; + if(begin == end) return; + } else if(available_bits % 32 == 0) { + for(int i = 0; i < available_bits / 32; ++i) { + boost::uint_least32_t word = boost::uint_least32_t(val) & 0xFFFFFFFFu; + int supress_warning = (bits >= 32); + BOOST_ASSERT(supress_warning == 1); + val >>= (32 * supress_warning); + *begin++ = word; + if(begin == end) return; + } + } else if(bits < 32 && available_bits < 32 - remaining_bits) { + saved_bits |= boost::uint_least32_t(val) << remaining_bits; + remaining_bits += bits; + } else { + if(bits < 32 || remaining_bits != 0) { + boost::uint_least32_t extra_bits = boost::uint_least32_t(val) & ((boost::uint_least32_t(1) << (32 - remaining_bits)) - 1); + val >>= 32 - remaining_bits; + *begin++ = saved_bits | (extra_bits << remaining_bits); + if(begin == end) return; + available_bits -= 32 - remaining_bits; + remaining_bits = 0; + } + if(bits >= 32) { + for(; available_bits >= 32; available_bits -= 32) { + boost::uint_least32_t word = boost::uint_least32_t(val) & 0xFFFFFFFFu; + int supress_warning = (bits >= 32); + BOOST_ASSERT(supress_warning == 1); + val >>= (32 * supress_warning); + *begin++ = word; + if(begin == end) return; + } + } + remaining_bits = available_bits; + saved_bits = static_cast<boost::uint_least32_t>(val); + } + } +} + +template<class Engine, class Iter> +void generate_impl(Engine& eng, Iter first, Iter last, boost::mpl::true_) +{ + return detail::generate_from_int(eng, first, last); +} + +template<class Engine, class Iter> +void generate_impl(Engine& eng, Iter first, Iter last, boost::mpl::false_) +{ + return detail::generate_from_real(eng, first, last); +} + +template<class Engine, class Iter> +void generate(Engine& eng, Iter first, Iter last) +{ + return detail::generate_impl(eng, first, last, boost::is_integral<typename Engine::result_type>()); +} + + + +template<class IntType, IntType m, class SeedSeq> +IntType seed_one_int(SeedSeq& seq) +{ + static const int log = ::boost::mpl::if_c<(m == 0), + ::boost::mpl::int_<(::std::numeric_limits<IntType>::digits)>, + ::boost::static_log2<m> >::type::value; + static const int k = + (log + ((~(static_cast<IntType>(2) << (log - 1)) & m)? 32 : 31)) / 32; + ::boost::uint_least32_t array[log / 32 + 4]; + seq.generate(&array[0], &array[0] + k + 3); + IntType s = 0; + for(int j = 0; j < k; ++j) { + IntType digit = const_mod<IntType, m>::apply(IntType(array[j+3])); + IntType mult = IntType(1) << 32*j; + s = const_mod<IntType, m>::mult_add(mult, digit, s); + } + return s; +} + +template<class IntType, IntType m, class Iter> +IntType get_one_int(Iter& first, Iter last) +{ + static const int log = ::boost::mpl::if_c<(m == 0), + ::boost::mpl::int_<(::std::numeric_limits<IntType>::digits)>, + ::boost::static_log2<m> >::type::value; + static const int k = + (log + ((~(static_cast<IntType>(2) << (log - 1)) & m)? 32 : 31)) / 32; + IntType s = 0; + for(int j = 0; j < k; ++j) { + if(first == last) { + throw ::std::invalid_argument("Not enough elements in call to seed."); + } + IntType digit = const_mod<IntType, m>::apply(IntType(*first++)); + IntType mult = IntType(1) << 32*j; + s = const_mod<IntType, m>::mult_add(mult, digit, s); + } + return s; +} + +// TODO: work in-place whenever possible +template<int w, std::size_t n, class SeedSeq, class UIntType> +void seed_array_int_impl(SeedSeq& seq, UIntType (&x)[n]) +{ + 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 += static_cast<UIntType>(storage[(w+31)/32*j + k]) << 32*k; + } + x[j] = val & ::boost::low_bits_mask_t<w>::sig_bits; + } +} + +template<int w, std::size_t n, class SeedSeq, class IntType> +inline void seed_array_int_impl(SeedSeq& seq, IntType (&x)[n], boost::mpl::true_) +{ + typedef typename boost::make_unsigned<IntType>::type unsigned_array[n]; + seed_array_int_impl<w>(seq, reinterpret_cast<unsigned_array&>(x)); +} + +template<int w, std::size_t n, class SeedSeq, class IntType> +inline void seed_array_int_impl(SeedSeq& seq, IntType (&x)[n], boost::mpl::false_) +{ + seed_array_int_impl<w>(seq, x); +} + +template<int w, std::size_t n, class SeedSeq, class IntType> +inline void seed_array_int(SeedSeq& seq, IntType (&x)[n]) +{ + seed_array_int_impl<w>(seq, x, boost::is_signed<IntType>()); +} + +template<int w, std::size_t n, class Iter, class UIntType> +void fill_array_int_impl(Iter& first, Iter last, UIntType (&x)[n]) +{ + for(std::size_t j = 0; j < n; j++) { + UIntType val = 0; + for(std::size_t k = 0; k < (w+31)/32; ++k) { + if(first == last) { + throw std::invalid_argument("Not enough elements in call to seed."); + } + val += static_cast<UIntType>(*first++) << 32*k; + } + x[j] = val & ::boost::low_bits_mask_t<w>::sig_bits; + } +} + +template<int w, std::size_t n, class Iter, class IntType> +inline void fill_array_int_impl(Iter& first, Iter last, IntType (&x)[n], boost::mpl::true_) +{ + typedef typename boost::make_unsigned<IntType>::type unsigned_array[n]; + fill_array_int_impl<w>(first, last, reinterpret_cast<unsigned_array&>(x)); +} + +template<int w, std::size_t n, class Iter, class IntType> +inline void fill_array_int_impl(Iter& first, Iter last, IntType (&x)[n], boost::mpl::false_) +{ + fill_array_int_impl<w>(first, last, x); +} + +template<int w, std::size_t n, class Iter, class IntType> +inline void fill_array_int(Iter& first, Iter last, IntType (&x)[n]) +{ + fill_array_int_impl<w>(first, last, x, boost::is_signed<IntType>()); +} + +template<int w, std::size_t n, class RealType> +void seed_array_real_impl(const boost::uint_least32_t* storage, RealType (&x)[n]) +{ + boost::uint_least32_t mask = ~((~boost::uint_least32_t(0)) << (w%32)); + RealType two32 = 4294967296.0; + const RealType divisor = RealType(1)/detail::pow2<RealType>(w); + unsigned int j; + for(j = 0; j < n; ++j) { + RealType val = RealType(0); + RealType mult = divisor; + for(int k = 0; k < w/32; ++k) { + val += *storage++ * mult; + mult *= two32; + } + if(mask != 0) { + val += (*storage++ & mask) * mult; + } + BOOST_ASSERT(val >= 0); + BOOST_ASSERT(val < 1); + x[j] = val; + } +} + +template<int w, std::size_t n, class SeedSeq, class RealType> +void seed_array_real(SeedSeq& seq, RealType (&x)[n]) +{ + using std::pow; + boost::uint_least32_t storage[((w+31)/32) * n]; + seq.generate(&storage[0], &storage[0] + ((w+31)/32) * n); + seed_array_real_impl<w>(storage, x); +} + +template<int w, std::size_t n, class Iter, class RealType> +void fill_array_real(Iter& first, Iter last, RealType (&x)[n]) +{ + boost::uint_least32_t mask = ~((~boost::uint_least32_t(0)) << (w%32)); + RealType two32 = 4294967296.0; + const RealType divisor = RealType(1)/detail::pow2<RealType>(w); + unsigned int j; + for(j = 0; j < n; ++j) { + RealType val = RealType(0); + RealType mult = divisor; + for(int k = 0; k < w/32; ++k, ++first) { + if(first == last) throw std::invalid_argument("Not enough elements in call to seed."); + val += *first * mult; + mult *= two32; + } + if(mask != 0) { + if(first == last) throw std::invalid_argument("Not enough elements in call to seed."); + val += (*first & mask) * mult; + ++first; + } + BOOST_ASSERT(val >= 0); + BOOST_ASSERT(val < 1); + x[j] = val; + } +} + +} +} +} + +#include <boost/random/detail/enable_warnings.hpp> + +#endif diff --git a/3rdParty/Boost/src/boost/random/detail/uniform_int_float.hpp b/3rdParty/Boost/src/boost/random/detail/uniform_int_float.hpp index 4607021..ef20915 100644 --- a/3rdParty/Boost/src/boost/random/detail/uniform_int_float.hpp +++ b/3rdParty/Boost/src/boost/random/detail/uniform_int_float.hpp @@ -1,85 +1,76 @@ /* boost random/detail/uniform_int_float.hpp header file * * Copyright Jens Maurer 2000-2001 + * Copyright Steven Watanabe 2011 * 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: uniform_int_float.hpp 52492 2009-04-19 14:55:57Z steven_watanabe $ + * $Id: uniform_int_float.hpp 72951 2011-07-07 04:57:37Z steven_watanabe $ * */ #ifndef BOOST_RANDOM_DETAIL_UNIFORM_INT_FLOAT_HPP #define BOOST_RANDOM_DETAIL_UNIFORM_INT_FLOAT_HPP +#include <boost/limits.hpp> #include <boost/config.hpp> +#include <boost/integer.hpp> #include <boost/random/detail/config.hpp> -#include <boost/random/uniform_01.hpp> +#include <boost/random/detail/generator_bits.hpp> +#include <boost/random/detail/disable_warnings.hpp> namespace boost { namespace random { namespace detail { -template<class UniformRandomNumberGenerator, class IntType = unsigned long> +template<class URNG> class uniform_int_float { public: - typedef UniformRandomNumberGenerator base_type; - typedef IntType result_type; + typedef URNG base_type; + typedef typename base_type::result_type base_result; - uniform_int_float(base_type rng, IntType min_arg = 0, IntType max_arg = 0xffffffff) - : _rng(rng), _min(min_arg), _max(max_arg) - { - init(); - } + typedef typename boost::uint_t< + (std::numeric_limits<boost::uintmax_t>::digits < + std::numeric_limits<base_result>::digits)? + std::numeric_limits<boost::uintmax_t>::digits : + std::numeric_limits<base_result>::digits + >::fast result_type; - result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return _min; } - result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return _max; } - base_type& base() { return _rng.base(); } - const base_type& base() const { return _rng.base(); } + uniform_int_float(base_type& rng) + : _rng(rng) {} - result_type operator()() - { - return static_cast<IntType>(_rng() * _range) + _min; - } + static result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () + { return 0; } + static result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () + { + std::size_t digits = std::numeric_limits<result_type>::digits; + if(detail::generator_bits<URNG>::value() < digits) { + digits = detail::generator_bits<URNG>::value(); + } + return (result_type(2) << (digits - 1)) - 1; + } + base_type& base() { return _rng; } + const base_type& base() const { return _rng; } -#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS - template<class CharT, class Traits> - friend std::basic_ostream<CharT,Traits>& - operator<<(std::basic_ostream<CharT,Traits>& os, const uniform_int_float& ud) - { - os << ud._min << " " << ud._max; - return os; - } - - template<class CharT, class Traits> - friend std::basic_istream<CharT,Traits>& - operator>>(std::basic_istream<CharT,Traits>& is, uniform_int_float& ud) - { - is >> std::ws >> ud._min >> std::ws >> ud._max; - ud.init(); - return is; - } -#endif + result_type operator()() + { + base_result range = static_cast<base_result>((max)())+1; + return static_cast<result_type>(_rng() * range); + } private: - void init() - { - _range = static_cast<base_result>(_max-_min)+1; - } - - typedef typename base_type::result_type base_result; - uniform_01<base_type> _rng; - result_type _min, _max; - base_result _range; + base_type& _rng; }; - } // namespace detail } // namespace random } // namespace boost +#include <boost/random/detail/enable_warnings.hpp> + #endif // BOOST_RANDOM_DETAIL_UNIFORM_INT_FLOAT_HPP diff --git a/3rdParty/Boost/src/boost/random/mersenne_twister.hpp b/3rdParty/Boost/src/boost/random/mersenne_twister.hpp index fa80aa6..78d0ab0 100644 --- a/3rdParty/Boost/src/boost/random/mersenne_twister.hpp +++ b/3rdParty/Boost/src/boost/random/mersenne_twister.hpp @@ -1,13 +1,14 @@ /* 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) * * See http://www.boost.org for most recent version including documentation. * - * $Id: mersenne_twister.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * $Id: mersenne_twister.hpp 71018 2011-04-05 21:27:52Z steven_watanabe $ * * Revision history * 2001-02-18 moved to individual header files @@ -16,25 +17,23 @@ #ifndef BOOST_RANDOM_MERSENNE_TWISTER_HPP #define BOOST_RANDOM_MERSENNE_TWISTER_HPP -#include <iostream> -#include <algorithm> // std::copy +#include <iosfwd> +#include <istream> #include <stdexcept> #include <boost/config.hpp> -#include <boost/limits.hpp> -#include <boost/static_assert.hpp> -#include <boost/integer_traits.hpp> #include <boost/cstdint.hpp> -#include <boost/random/linear_congruential.hpp> -#include <boost/detail/workaround.hpp> +#include <boost/integer/integer_mask.hpp> #include <boost/random/detail/config.hpp> #include <boost/random/detail/ptr_helper.hpp> #include <boost/random/detail/seed.hpp> +#include <boost/random/detail/seed_impl.hpp> +#include <boost/random/detail/generator_seed_seq.hpp> 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,277 +60,406 @@ namespace random { * its state array. For example, \mt11213b requires about 1408 bytes and * \mt19937 requires about 2496 bytes. */ -template<class UIntType, int w, int n, int m, int r, UIntType a, int u, - int s, UIntType b, int t, UIntType c, int l, UIntType val> -class mersenne_twister +template<class UIntType, + std::size_t w, std::size_t n, std::size_t m, std::size_t r, + UIntType a, std::size_t u, UIntType d, std::size_t s, + UIntType b, std::size_t t, + UIntType c, std::size_t l, UIntType f> +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(); } - - /** - * Constructs a @c mersenne_twister and calls @c seed(value). - */ - BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(mersenne_twister, UIntType, value) - { seed(value); } - template<class It> 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); } - - // compiler-generated copy ctor and assignment operator are fine - - /** Calls @c seed(result_type(5489)). */ - void seed() { seed(UIntType(5489)); } - - /** - * 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 2<sup>w</sup> - * 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; + // 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_engine and calls @c seed(). + */ + mersenne_twister_engine() { seed(); } + + /** + * Constructs a @c mersenne_twister_engine and calls @c seed(value). + */ + BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(mersenne_twister_engine, + UIntType, value) + { seed(value); } + template<class It> mersenne_twister_engine(It& first, It last) + { seed(first,last); } + + /** + * 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); } + + // 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 2<sup>w</sup> + * 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) + { + detail::seed_array_int<w>(seq, x); + i = n; + + // fix up the state if it's all zeroes. + if((x[0] & (~static_cast<UIntType>(0) << r)) == 0) { + for(std::size_t j = 1; i < n; ++j) { + if(x[j] != 0) return; + } + x[0] = static_cast<UIntType>(1) << (w-1); + } } - } - - /** - * 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<result_type>::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; - } - - template<class It> - 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"); - } - - 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; - } - - result_type operator()(); - static bool validation(result_type v) { return val == v; } -#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE + /** Sets the state of the generator using values from an iterator range. */ + template<class It> + void seed(It& first, It last) + { + detail::fill_array_int<w>(first, last, x); + i = n; + + // fix up the state if it's all zeroes. + if((x[0] & (~static_cast<UIntType>(0) << r)) == 0) { + for(std::size_t j = 1; i < n; ++j) { + if(x[j] != 0) return; + } + x[0] = static_cast<UIntType>(1) << (w-1); + } + } + + /** 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 () + { return boost::low_bits_mask_t<w>::sig_bits; } + + /** Produces the next value of the generator. */ + result_type operator()(); + + /** Fills a range with random values */ + template<class Iter> + void generate(Iter first, Iter last) + { detail::generate_from_int(*this, first, last); } + + /** + * 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::uintmax_t z) + { + for(boost::uintmax_t j = 0; j < z; ++j) { + (*this)(); + } + } #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS - template<class CharT, class Traits> - friend std::basic_ostream<CharT,Traits>& - operator<<(std::basic_ostream<CharT,Traits>& os, const mersenne_twister& mt) - { - for(int j = 0; j < mt.state_size; ++j) - os << mt.compute(j) << " "; - return os; - } - - template<class CharT, class Traits> - friend std::basic_istream<CharT,Traits>& - operator>>(std::basic_istream<CharT,Traits>& 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<class CharT, class Traits> + friend std::basic_ostream<CharT,Traits>& + operator<<(std::basic_ostream<CharT,Traits>& os, + const mersenne_twister_engine& mt) + { + mt.print(os); + return os; + } + + /** Reads a mersenne_twister_engine from a @c std::istream */ + template<class CharT, class Traits> + friend std::basic_istream<CharT,Traits>& + operator>>(std::basic_istream<CharT,Traits>& 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) + { + if(x.i < y.i) return x.equal_imp(y); + else return y.equal_imp(x); + } + + /** + * 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 - - // 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; + /// \cond show_private + + void twist(); + + /** + * Does the work of operator==. This is in a member function + * for portability. Some compilers, such as msvc 7.1 and + * Sun CC 5.10 can't access template parameters or static + * members of the class from inline friend functions. + * + * requires i <= other.i + */ + bool equal_imp(const mersenne_twister_engine& other) const + { + UIntType back[n]; + std::size_t offset = other.i - i; + for(std::size_t j = 0; j + offset < n; ++j) + if(x[j] != other.x[j+offset]) + return false; + rewind(&back[n-1], offset); + for(std::size_t j = 0; j < offset; ++j) + if(back[j + n - offset] != other.x[j]) + return false; + return true; + } + + /** + * Does the work of operator<<. This is in a member function + * for portability. + */ + template<class CharT, class Traits> + void print(std::basic_ostream<CharT, Traits>& os) const + { + UIntType data[n]; + for(std::size_t j = 0; j < i; ++j) { + data[j + n - i] = x[j]; + } + if(i != n) { + rewind(&data[n - i - 1], n - i); + } + os << data[0]; + for(std::size_t j = 1; j < n; ++j) { + os << ' ' << data[j]; + } + } + + /** + * Copies z elements of the state preceding x[0] into + * the array whose last element is last. + */ + void rewind(UIntType* last, std::size_t z) const + { + const UIntType upper_mask = (~static_cast<UIntType>(0)) << r; + const UIntType lower_mask = ~upper_mask; + UIntType y0 = x[m-1] ^ x[n-1]; + if(y0 & (static_cast<UIntType>(1) << (w-1))) { + y0 = ((y0 ^ a) << 1) | 1; + } else { + y0 = y0 << 1; + } + for(std::size_t sz = 0; sz < z; ++sz) { + UIntType y1 = + rewind_find(last, sz, m-1) ^ rewind_find(last, sz, n-1); + if(y1 & (static_cast<UIntType>(1) << (w-1))) { + y1 = ((y1 ^ a) << 1) | 1; + } else { + y1 = y1 << 1; + } + *(last - sz) = (y0 & upper_mask) | (y1 & lower_mask); + y0 = y1; + } + } + + /** + * Given a pointer to the last element of the rewind array, + * and the current size of the rewind array, finds an element + * relative to the next available slot in the rewind array. + */ + UIntType + rewind_find(UIntType* last, std::size_t size, std::size_t j) const + { + std::size_t index = (j + n - size + n - 1) % n; + if(index < n - size) { + return x[index]; + } else { + return *(last - (n - 1 - index)); + } + } + + /// \endcond + + // state representation: next output is o(x(i)) + // x[0] ... x[k] x[k+1] ... x[n-1] represents + // x(i-k) ... x(i) x(i+1) ... x(i-k+n-1) + + UIntType x[n]; + std::size_t i; }; +/// \cond show_private + #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION // A definition is required even for integral static constants -template<class UIntType, int w, int n, int m, int r, UIntType a, int u, - int s, UIntType b, int t, UIntType c, int l, UIntType val> -const bool mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::has_fixed_range; -template<class UIntType, int w, int n, int m, int r, UIntType a, int u, - int s, UIntType b, int t, UIntType c, int l, UIntType val> -const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::state_size; -template<class UIntType, int w, int n, int m, int r, UIntType a, int u, - int s, UIntType b, int t, UIntType c, int l, UIntType val> -const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::shift_size; -template<class UIntType, int w, int n, int m, int r, UIntType a, int u, - int s, UIntType b, int t, UIntType c, int l, UIntType val> -const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::mask_bits; -template<class UIntType, int w, int n, int m, int r, UIntType a, int u, - int s, UIntType b, int t, UIntType c, int l, UIntType val> -const UIntType mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::parameter_a; -template<class UIntType, int w, int n, int m, int r, UIntType a, int u, - int s, UIntType b, int t, UIntType c, int l, UIntType val> -const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_u; -template<class UIntType, int w, int n, int m, int r, UIntType a, int u, - int s, UIntType b, int t, UIntType c, int l, UIntType val> -const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_s; -template<class UIntType, int w, int n, int m, int r, UIntType a, int u, - int s, UIntType b, int t, UIntType c, int l, UIntType val> -const UIntType mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_b; -template<class UIntType, int w, int n, int m, int r, UIntType a, int u, - int s, UIntType b, int t, UIntType c, int l, UIntType val> -const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_t; -template<class UIntType, int w, int n, int m, int r, UIntType a, int u, - int s, UIntType b, int t, UIntType c, int l, UIntType val> -const UIntType mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_c; -template<class UIntType, int w, int n, int m, int r, UIntType a, int u, - int s, UIntType b, int t, UIntType c, int l, UIntType val> -const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_l; +#define BOOST_RANDOM_MT_DEFINE_CONSTANT(type, name) \ +template<class UIntType, std::size_t w, std::size_t n, std::size_t m, \ + std::size_t r, UIntType a, std::size_t u, UIntType d, std::size_t s, \ + UIntType b, std::size_t t, UIntType c, std::size_t l, UIntType f> \ +const type mersenne_twister_engine<UIntType,w,n,m,r,a,u,d,s,b,t,c,l,f>::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<class UIntType, int w, int n, int m, int r, UIntType a, int u, - int s, UIntType b, int t, UIntType c, int l, UIntType val> -void mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::twist(int block) +template<class UIntType, + std::size_t w, std::size_t n, std::size_t m, std::size_t r, + UIntType a, std::size_t u, UIntType d, std::size_t s, + UIntType b, std::size_t t, + UIntType c, std::size_t l, UIntType f> +void +mersenne_twister_engine<UIntType,w,n,m,r,a,u,d,s,b,t,c,l,f>::twist() { - const UIntType upper_mask = (~0u) << r; - const UIntType lower_mask = ~upper_mask; + const UIntType upper_mask = (~static_cast<UIntType>(0)) << r; + const UIntType lower_mask = ~upper_mask; + + const std::size_t unroll_factor = 6; + const std::size_t unroll_extra1 = (n-m) % unroll_factor; + const std::size_t unroll_extra2 = (m-1) % unroll_factor; - 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); - } + for(std::size_t j = 0; j < n-m-unroll_extra1; j++) { + UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask); + x[j] = x[j+m] ^ (y >> 1) ^ ((x[j+1]&1) * a); + } } - - 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-unroll_extra1; j < n-m; j++) { + UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask); + x[j] = x[j+m] ^ (y >> 1) ^ ((x[j+1]&1) * a); + } + } + { + for(std::size_t j = n-m; j < n-1-unroll_extra2; j++) { + UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask); + x[j] = x[j-(n-m)] ^ (y >> 1) ^ ((x[j+1]&1) * a); + } + } + { + for(std::size_t j = n-1-unroll_extra2; j < n-1; j++) { + UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask); + x[j] = x[j-(n-m)] ^ (y >> 1) ^ ((x[j+1]&1) * a); + } } // 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); + UIntType y = (x[n-1] & upper_mask) | (x[0] & lower_mask); + x[n-1] = x[m-1] ^ (y >> 1) ^ ((x[0]&1) * a); i = 0; - } } /// \endcond -template<class UIntType, int w, int n, int m, int r, UIntType a, int u, - int s, UIntType b, int t, UIntType c, int l, UIntType val> -inline typename mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::result_type -mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::operator()() +template<class UIntType, + std::size_t w, std::size_t n, std::size_t m, std::size_t r, + UIntType a, std::size_t u, UIntType d, std::size_t s, + UIntType b, std::size_t t, + UIntType c, std::size_t l, UIntType f> +inline typename +mersenne_twister_engine<UIntType,w,n,m,r,a,u,d,s,b,t,c,l,f>::result_type +mersenne_twister_engine<UIntType,w,n,m,r,a,u,d,s,b,t,c,l,f>::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(); + // Step 4 + UIntType z = x[i]; + ++i; + z ^= ((z >> u) & d); + z ^= ((z << s) & b); + z ^= ((z << t) & c); + z ^= (z >> l); + return z; } -} // namespace random - /** * The specializations \mt11213b and \mt19937 are from * @@ -343,8 +471,8 @@ mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::operator()() * Generation, Vol. 8, No. 1, January 1998, pp. 3-30. * @endblockquote */ -typedef random::mersenne_twister<uint32_t,32,351,175,19,0xccab8ee7,11, - 7,0x31b6ab00,15,0xffe50000,17, 0xa37d3c92> mt11213b; +typedef mersenne_twister_engine<uint32_t,32,351,175,19,0xccab8ee7, + 11,0xffffffff,7,0x31b6ab00,15,0xffe50000,17,1812433253> mt11213b; /** * The specializations \mt11213b and \mt19937 are from @@ -357,11 +485,61 @@ typedef random::mersenne_twister<uint32_t,32,351,175,19,0xccab8ee7,11, * Generation, Vol. 8, No. 1, January 1998, pp. 3-30. * @endblockquote */ -typedef random::mersenne_twister<uint32_t,32,624,397,31,0x9908b0df,11, - 7,0x9d2c5680,15,0xefc60000,18, 3346425566U> mt19937; +typedef mersenne_twister_engine<uint32_t,32,624,397,31,0x9908b0df, + 11,0xffffffff,7,0x9d2c5680,15,0xefc60000,18,1812433253> mt19937; + +#if !defined(BOOST_NO_INT64_T) && !defined(BOOST_NO_INTEGRAL_INT64_T) +typedef mersenne_twister_engine<uint64_t,64,312,156,31, + UINT64_C(0xb5026f5aa96619e9),29,UINT64_C(0x5555555555555555),17, + UINT64_C(0x71d67fffeda60000),37,UINT64_C(0xfff7eee000000000),43, + UINT64_C(6364136223846793005)> mt19937_64; +#endif + +/// \cond show_deprecated + +template<class UIntType, + int w, int n, int m, int r, + UIntType a, int u, std::size_t s, + UIntType b, int t, + UIntType c, int l, UIntType v> +class mersenne_twister : + public mersenne_twister_engine<UIntType, + w, n, m, r, a, u, ~(UIntType)0, s, b, t, c, l, 1812433253> +{ + typedef mersenne_twister_engine<UIntType, + w, n, m, r, a, u, ~(UIntType)0, s, b, t, c, l, 1812433253> base_type; +public: + mersenne_twister() {} + BOOST_RANDOM_DETAIL_GENERATOR_CONSTRUCTOR(mersenne_twister, Gen, gen) + { seed(gen); } + BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(mersenne_twister, UIntType, val) + { seed(val); } + template<class It> + mersenne_twister(It& first, It last) : base_type(first, last) {} + void seed() { base_type::seed(); } + BOOST_RANDOM_DETAIL_GENERATOR_SEED(mersenne_twister, Gen, gen) + { + detail::generator_seed_seq<Gen> seq(gen); + base_type::seed(seq); + } + BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(mersenne_twister, UIntType, val) + { base_type::seed(val); } + template<class It> + void seed(It& first, It last) { base_type::seed(first, last); } +}; + +/// \endcond + +} // namespace random + +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/3rdParty/Boost/src/boost/random/uniform_int.hpp b/3rdParty/Boost/src/boost/random/uniform_int.hpp index 426a9e1..7ae3b92 100644 --- a/3rdParty/Boost/src/boost/random/uniform_int.hpp +++ b/3rdParty/Boost/src/boost/random/uniform_int.hpp @@ -7,7 +7,7 @@ * * See http://www.boost.org for most recent version including documentation. * - * $Id: uniform_int.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * $Id: uniform_int.hpp 71018 2011-04-05 21:27:52Z steven_watanabe $ * * Revision history * 2001-04-08 added min<max assertion (N. Becker) @@ -17,15 +17,8 @@ #ifndef BOOST_RANDOM_UNIFORM_INT_HPP #define BOOST_RANDOM_UNIFORM_INT_HPP -#include <cassert> -#include <iostream> -#include <boost/config.hpp> -#include <boost/limits.hpp> -#include <boost/static_assert.hpp> -#include <boost/detail/workaround.hpp> -#include <boost/random/detail/config.hpp> -#include <boost/random/detail/signed_unsigned_tools.hpp> -#include <boost/type_traits/make_unsigned.hpp> +#include <boost/assert.hpp> +#include <boost/random/uniform_int_distribution.hpp> namespace boost { @@ -35,264 +28,70 @@ namespace boost { * distributed in the set of integer numbers {min, min+1, min+2, ..., max}. * * The template parameter IntType shall denote an integer-like value type. + * + * This class is deprecated. Please use @c uniform_int_distribution in + * new code. */ template<class IntType = int> -class uniform_int +class uniform_int : public random::uniform_int_distribution<IntType> { + typedef random::uniform_int_distribution<IntType> base_type; public: - typedef IntType input_type; - typedef IntType result_type; - - /// \cond hide_private_members - typedef typename make_unsigned<result_type>::type range_type; - /// \endcond - - /** - * Constructs a uniform_int object. @c min and @c max are - * the parameters of the distribution. - * - * Requires: min <= max - */ - explicit uniform_int(IntType min_arg = 0, IntType max_arg = 9) - : _min(min_arg), _max(max_arg) - { -#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS - // MSVC fails BOOST_STATIC_ASSERT with std::numeric_limits at class scope - BOOST_STATIC_ASSERT(std::numeric_limits<IntType>::is_integer); -#endif - assert(min_arg <= max_arg); - init(); - } - - /** - * Returns: The "min" parameter of the distribution - */ - result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return _min; } - /** - * Returns: The "max" parameter of the distribution - */ - result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return _max; } - void reset() { } - - // can't have member function templates out-of-line due to MSVC bugs - template<class Engine> - result_type operator()(Engine& eng) - { - return generate(eng, _min, _max, _range); - } - - template<class Engine> - result_type operator()(Engine& eng, result_type n) - { - assert(n > 0); - - if (n == 1) - { - return 0; - } - - return generate(eng, 0, n - 1, n - 1); - } - -#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS - template<class CharT, class Traits> - friend std::basic_ostream<CharT,Traits>& - operator<<(std::basic_ostream<CharT,Traits>& os, const uniform_int& ud) - { - os << ud._min << " " << ud._max; - return os; - } - - template<class CharT, class Traits> - friend std::basic_istream<CharT,Traits>& - operator>>(std::basic_istream<CharT,Traits>& is, uniform_int& ud) - { - is >> std::ws >> ud._min >> std::ws >> ud._max; - ud.init(); - return is; - } -#endif - -private: - -#ifdef BOOST_MSVC -#pragma warning(push) -// disable division by zero warning, since we can't -// actually divide by zero. -#pragma warning(disable:4723) -#endif - /// \cond hide_private_members - template<class Engine> - static result_type generate(Engine& eng, result_type min_value, result_type /*max_value*/, range_type range) - { - typedef typename Engine::result_type base_result; - // ranges are always unsigned - typedef typename make_unsigned<base_result>::type base_unsigned; - const base_result bmin = (eng.min)(); - const base_unsigned brange = - random::detail::subtract<base_result>()((eng.max)(), (eng.min)()); - - if(range == 0) { - return min_value; - } else if(brange == range) { - // this will probably never happen in real life - // basically nothing to do; just take care we don't overflow / underflow - base_unsigned v = random::detail::subtract<base_result>()(eng(), bmin); - return random::detail::add<base_unsigned, result_type>()(v, min_value); - } else if(brange < range) { - // use rejection method to handle things like 0..3 --> 0..4 - for(;;) { - // concatenate several invocations of the base RNG - // take extra care to avoid overflows - - // limit == floor((range+1)/(brange+1)) - // Therefore limit*(brange+1) <= range+1 - range_type limit; - if(range == (std::numeric_limits<range_type>::max)()) { - limit = range/(range_type(brange)+1); - if(range % (range_type(brange)+1) == range_type(brange)) - ++limit; - } else { - limit = (range+1)/(range_type(brange)+1); - } - - // We consider "result" as expressed to base (brange+1): - // For every power of (brange+1), we determine a random factor - range_type result = range_type(0); - range_type mult = range_type(1); - - // loop invariants: - // result < mult - // mult <= range - while(mult <= limit) { - // Postcondition: result <= range, thus no overflow - // - // limit*(brange+1)<=range+1 def. of limit (1) - // eng()-bmin<=brange eng() post. (2) - // and mult<=limit. loop condition (3) - // Therefore mult*(eng()-bmin+1)<=range+1 by (1),(2),(3) (4) - // Therefore mult*(eng()-bmin)+mult<=range+1 rearranging (4) (5) - // result<mult loop invariant (6) - // Therefore result+mult*(eng()-bmin)<range+1 by (5), (6) (7) - // - // Postcondition: result < mult*(brange+1) - // - // result<mult loop invariant (1) - // eng()-bmin<=brange eng() post. (2) - // Therefore result+mult*(eng()-bmin) < - // mult+mult*(eng()-bmin) by (1) (3) - // Therefore result+(eng()-bmin)*mult < - // mult+mult*brange by (2), (3) (4) - // Therefore result+(eng()-bmin)*mult < - // mult*(brange+1) by (4) - result += static_cast<range_type>(random::detail::subtract<base_result>()(eng(), bmin) * mult); - - // equivalent to (mult * (brange+1)) == range+1, but avoids overflow. - if(mult * range_type(brange) == range - mult + 1) { - // The destination range is an integer power of - // the generator's range. - return(result); - } - - // Postcondition: mult <= range - // - // limit*(brange+1)<=range+1 def. of limit (1) - // mult<=limit loop condition (2) - // Therefore mult*(brange+1)<=range+1 by (1), (2) (3) - // mult*(brange+1)!=range+1 preceding if (4) - // Therefore mult*(brange+1)<range+1 by (3), (4) (5) - // - // Postcondition: result < mult - // - // See the second postcondition on the change to result. - mult *= range_type(brange)+range_type(1); - } - // loop postcondition: range/mult < brange+1 - // - // mult > limit loop condition (1) - // Suppose range/mult >= brange+1 Assumption (2) - // range >= mult*(brange+1) by (2) (3) - // range+1 > mult*(brange+1) by (3) (4) - // range+1 > (limit+1)*(brange+1) by (1), (4) (5) - // (range+1)/(brange+1) > limit+1 by (5) (6) - // limit < floor((range+1)/(brange+1)) by (6) (7) - // limit==floor((range+1)/(brange+1)) def. of limit (8) - // not (2) reductio (9) - // - // loop postcondition: (range/mult)*mult+(mult-1) >= range - // - // (range/mult)*mult + range%mult == range identity (1) - // range%mult < mult def. of % (2) - // (range/mult)*mult+mult > range by (1), (2) (3) - // (range/mult)*mult+(mult-1) >= range by (3) (4) - // - // Note that the maximum value of result at this point is (mult-1), - // so after this final step, we generate numbers that can be - // at least as large as range. We have to really careful to avoid - // overflow in this final addition and in the rejection. Anything - // that overflows is larger than range and can thus be rejected. - - // range/mult < brange+1 -> no endless loop - range_type result_increment = uniform_int<range_type>(0, range/mult)(eng); - if((std::numeric_limits<range_type>::max)() / mult < result_increment) { - // The multiplcation would overflow. Reject immediately. - continue; - } - result_increment *= mult; - // unsigned integers are guaranteed to wrap on overflow. - result += result_increment; - if(result < result_increment) { - // The addition overflowed. Reject. - continue; - } - if(result > range) { - // Too big. Reject. - continue; - } - return random::detail::add<range_type, result_type>()(result, min_value); - } - } else { // brange > range - base_unsigned bucket_size; - // it's safe to add 1 to range, as long as we cast it first, - // because we know that it is less than brange. However, - // we do need to be careful not to cause overflow by adding 1 - // to brange. - if(brange == (std::numeric_limits<base_unsigned>::max)()) { - bucket_size = brange / (static_cast<base_unsigned>(range)+1); - if(brange % (static_cast<base_unsigned>(range)+1) == static_cast<base_unsigned>(range)) { - ++bucket_size; - } - } else { - bucket_size = (brange+1) / (static_cast<base_unsigned>(range)+1); - } - for(;;) { - base_unsigned result = - random::detail::subtract<base_result>()(eng(), bmin); - result /= bucket_size; - // result and range are non-negative, and result is possibly larger - // than range, so the cast is safe - if(result <= static_cast<base_unsigned>(range)) - return random::detail::add<base_unsigned, result_type>()(result, min_value); - } + class param_type : public base_type::param_type + { + public: + typedef uniform_int distribution_type; + /** + * Constructs the parameters of a uniform_int distribution. + * + * Requires: min <= max + */ + explicit param_type(IntType min_arg = 0, IntType max_arg = 9) + : base_type::param_type(min_arg, max_arg) + {} + }; + + /** + * Constructs a uniform_int object. @c min and @c max are + * the parameters of the distribution. + * + * Requires: min <= max + */ + explicit uniform_int(IntType min_arg = 0, IntType max_arg = 9) + : base_type(min_arg, max_arg) + {} + + /** Constructs a uniform_int distribution from its parameters. */ + explicit uniform_int(const param_type& parm) + : base_type(parm) + {} + + /** Returns the parameters of the distribution */ + param_type param() const { return param_type(this->a(), this->b()); } + /** Sets the parameters of the distribution. */ + void param(const param_type& parm) { this->base_type::param(parm); } + + // Codergear seems to have trouble with a using declaration here + + template<class Engine> + IntType operator()(Engine& eng) const + { + return static_cast<const base_type&>(*this)(eng); } - } - -#ifdef BOOST_MSVC -#pragma warning(pop) -#endif - - void init() - { - _range = random::detail::subtract<result_type>()(_max, _min); - } - /// \endcond + template<class Engine> + IntType operator()(Engine& eng, const param_type& parm) const + { + return static_cast<const base_type&>(*this)(eng, parm); + } - // The result_type may be signed or unsigned, but the _range is always - // unsigned. - result_type _min, _max; - range_type _range; + template<class Engine> + IntType operator()(Engine& eng, IntType n) const + { + BOOST_ASSERT(n > 0); + return static_cast<const base_type&>(*this)(eng, param_type(0, n - 1)); + } }; } // namespace boost diff --git a/3rdParty/Boost/src/boost/random/uniform_int_distribution.hpp b/3rdParty/Boost/src/boost/random/uniform_int_distribution.hpp new file mode 100644 index 0000000..0612028 --- /dev/null +++ b/3rdParty/Boost/src/boost/random/uniform_int_distribution.hpp @@ -0,0 +1,400 @@ +/* boost random/uniform_int_distribution.hpp header file + * + * Copyright Jens Maurer 2000-2001 + * Copyright Steven Watanabe 2011 + * 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: uniform_int_distribution.hpp 71018 2011-04-05 21:27:52Z steven_watanabe $ + * + * Revision history + * 2001-04-08 added min<max assertion (N. Becker) + * 2001-02-18 moved to individual header files + */ + +#ifndef BOOST_RANDOM_UNIFORM_INT_DISTRIBUTION_HPP +#define BOOST_RANDOM_UNIFORM_INT_DISTRIBUTION_HPP + +#include <iosfwd> +#include <ios> +#include <istream> +#include <boost/config.hpp> +#include <boost/limits.hpp> +#include <boost/assert.hpp> +#include <boost/random/detail/config.hpp> +#include <boost/random/detail/operators.hpp> +#include <boost/random/detail/uniform_int_float.hpp> +#include <boost/random/detail/signed_unsigned_tools.hpp> +#include <boost/type_traits/make_unsigned.hpp> +#include <boost/type_traits/is_integral.hpp> + +namespace boost { +namespace random { +namespace detail { + + +#ifdef BOOST_MSVC +#pragma warning(push) +// disable division by zero warning, since we can't +// actually divide by zero. +#pragma warning(disable:4723) +#endif + +template<class Engine, class T> +T generate_uniform_int( + Engine& eng, T min_value, T max_value, + boost::mpl::true_ /** is_integral<Engine::result_type> */) +{ + typedef T result_type; + typedef typename make_unsigned<T>::type range_type; + typedef typename Engine::result_type base_result; + // ranges are always unsigned + typedef typename make_unsigned<base_result>::type base_unsigned; + const range_type range = random::detail::subtract<result_type>()(max_value, min_value); + const base_result bmin = (eng.min)(); + const base_unsigned brange = + random::detail::subtract<base_result>()((eng.max)(), (eng.min)()); + + if(range == 0) { + return min_value; + } else if(brange == range) { + // this will probably never happen in real life + // basically nothing to do; just take care we don't overflow / underflow + base_unsigned v = random::detail::subtract<base_result>()(eng(), bmin); + return random::detail::add<base_unsigned, result_type>()(v, min_value); + } else if(brange < range) { + // use rejection method to handle things like 0..3 --> 0..4 + for(;;) { + // concatenate several invocations of the base RNG + // take extra care to avoid overflows + + // limit == floor((range+1)/(brange+1)) + // Therefore limit*(brange+1) <= range+1 + range_type limit; + if(range == (std::numeric_limits<range_type>::max)()) { + limit = range/(range_type(brange)+1); + if(range % (range_type(brange)+1) == range_type(brange)) + ++limit; + } else { + limit = (range+1)/(range_type(brange)+1); + } + + // We consider "result" as expressed to base (brange+1): + // For every power of (brange+1), we determine a random factor + range_type result = range_type(0); + range_type mult = range_type(1); + + // loop invariants: + // result < mult + // mult <= range + while(mult <= limit) { + // Postcondition: result <= range, thus no overflow + // + // limit*(brange+1)<=range+1 def. of limit (1) + // eng()-bmin<=brange eng() post. (2) + // and mult<=limit. loop condition (3) + // Therefore mult*(eng()-bmin+1)<=range+1 by (1),(2),(3) (4) + // Therefore mult*(eng()-bmin)+mult<=range+1 rearranging (4) (5) + // result<mult loop invariant (6) + // Therefore result+mult*(eng()-bmin)<range+1 by (5), (6) (7) + // + // Postcondition: result < mult*(brange+1) + // + // result<mult loop invariant (1) + // eng()-bmin<=brange eng() post. (2) + // Therefore result+mult*(eng()-bmin) < + // mult+mult*(eng()-bmin) by (1) (3) + // Therefore result+(eng()-bmin)*mult < + // mult+mult*brange by (2), (3) (4) + // Therefore result+(eng()-bmin)*mult < + // mult*(brange+1) by (4) + result += static_cast<range_type>(random::detail::subtract<base_result>()(eng(), bmin) * mult); + + // equivalent to (mult * (brange+1)) == range+1, but avoids overflow. + if(mult * range_type(brange) == range - mult + 1) { + // The destination range is an integer power of + // the generator's range. + return(result); + } + + // Postcondition: mult <= range + // + // limit*(brange+1)<=range+1 def. of limit (1) + // mult<=limit loop condition (2) + // Therefore mult*(brange+1)<=range+1 by (1), (2) (3) + // mult*(brange+1)!=range+1 preceding if (4) + // Therefore mult*(brange+1)<range+1 by (3), (4) (5) + // + // Postcondition: result < mult + // + // See the second postcondition on the change to result. + mult *= range_type(brange)+range_type(1); + } + // loop postcondition: range/mult < brange+1 + // + // mult > limit loop condition (1) + // Suppose range/mult >= brange+1 Assumption (2) + // range >= mult*(brange+1) by (2) (3) + // range+1 > mult*(brange+1) by (3) (4) + // range+1 > (limit+1)*(brange+1) by (1), (4) (5) + // (range+1)/(brange+1) > limit+1 by (5) (6) + // limit < floor((range+1)/(brange+1)) by (6) (7) + // limit==floor((range+1)/(brange+1)) def. of limit (8) + // not (2) reductio (9) + // + // loop postcondition: (range/mult)*mult+(mult-1) >= range + // + // (range/mult)*mult + range%mult == range identity (1) + // range%mult < mult def. of % (2) + // (range/mult)*mult+mult > range by (1), (2) (3) + // (range/mult)*mult+(mult-1) >= range by (3) (4) + // + // Note that the maximum value of result at this point is (mult-1), + // so after this final step, we generate numbers that can be + // at least as large as range. We have to really careful to avoid + // overflow in this final addition and in the rejection. Anything + // that overflows is larger than range and can thus be rejected. + + // range/mult < brange+1 -> no endless loop + range_type result_increment = + generate_uniform_int( + eng, + static_cast<range_type>(0), + static_cast<range_type>(range/mult), + boost::mpl::true_()); + if((std::numeric_limits<range_type>::max)() / mult < result_increment) { + // The multiplcation would overflow. Reject immediately. + continue; + } + result_increment *= mult; + // unsigned integers are guaranteed to wrap on overflow. + result += result_increment; + if(result < result_increment) { + // The addition overflowed. Reject. + continue; + } + if(result > range) { + // Too big. Reject. + continue; + } + return random::detail::add<range_type, result_type>()(result, min_value); + } + } else { // brange > range + base_unsigned bucket_size; + // it's safe to add 1 to range, as long as we cast it first, + // because we know that it is less than brange. However, + // we do need to be careful not to cause overflow by adding 1 + // to brange. + if(brange == (std::numeric_limits<base_unsigned>::max)()) { + bucket_size = brange / (static_cast<base_unsigned>(range)+1); + if(brange % (static_cast<base_unsigned>(range)+1) == static_cast<base_unsigned>(range)) { + ++bucket_size; + } + } else { + bucket_size = (brange+1) / (static_cast<base_unsigned>(range)+1); + } + for(;;) { + base_unsigned result = + random::detail::subtract<base_result>()(eng(), bmin); + result /= bucket_size; + // result and range are non-negative, and result is possibly larger + // than range, so the cast is safe + if(result <= static_cast<base_unsigned>(range)) + return random::detail::add<base_unsigned, result_type>()(result, min_value); + } + } +} + +#ifdef BOOST_MSVC +#pragma warning(pop) +#endif + +template<class Engine, class T> +inline T generate_uniform_int( + Engine& eng, T min_value, T max_value, + boost::mpl::false_ /** is_integral<Engine::result_type> */) +{ + uniform_int_float<Engine> wrapper(eng); + return generate_uniform_int(wrapper, min_value, max_value, boost::mpl::true_()); +} + +template<class Engine, class T> +inline T generate_uniform_int(Engine& eng, T min_value, T max_value) +{ + typedef typename Engine::result_type base_result; + return generate_uniform_int(eng, min_value, max_value, + boost::is_integral<base_result>()); +} + +} + +/** + * The class template uniform_int_distribution models a \random_distribution. + * On each invocation, it returns a random integer value uniformly + * distributed in the set of integers {min, min+1, min+2, ..., max}. + * + * The template parameter IntType shall denote an integer-like value type. + */ +template<class IntType = int> +class uniform_int_distribution +{ +public: + typedef IntType input_type; + typedef IntType result_type; + + class param_type + { + public: + + typedef uniform_int_distribution distribution_type; + + /** + * Constructs the parameters of a uniform_int_distribution. + * + * Requires min <= max + */ + explicit param_type( + IntType min_arg = 0, + IntType max_arg = (std::numeric_limits<IntType>::max)()) + : _min(min_arg), _max(max_arg) + { + BOOST_ASSERT(_min <= _max); + } + + /** Returns the minimum value of the distribution. */ + IntType a() const { return _min; } + /** Returns the maximum value of the distribution. */ + IntType b() const { return _max; } + + /** Writes the parameters to a @c std::ostream. */ + BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm) + { + os << parm._min << " " << parm._max; + return os; + } + + /** Reads the parameters from a @c std::istream. */ + BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm) + { + IntType min_in, max_in; + if(is >> min_in >> std::ws >> max_in) { + if(min_in <= max_in) { + parm._min = min_in; + parm._max = max_in; + } else { + is.setstate(std::ios_base::failbit); + } + } + return is; + } + + /** Returns true if the two sets of parameters are equal. */ + BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs) + { return lhs._min == rhs._min && lhs._max == rhs._max; } + + /** Returns true if the two sets of parameters are different. */ + BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type) + + private: + + IntType _min; + IntType _max; + }; + + /** + * Constructs a uniform_int_distribution. @c min and @c max are + * the parameters of the distribution. + * + * Requires: min <= max + */ + explicit uniform_int_distribution( + IntType min_arg = 0, + IntType max_arg = (std::numeric_limits<IntType>::max)()) + : _min(min_arg), _max(max_arg) + { + BOOST_ASSERT(min_arg <= max_arg); + } + /** Constructs a uniform_int_distribution from its parameters. */ + explicit uniform_int_distribution(const param_type& parm) + : _min(parm.a()), _max(parm.b()) {} + + /** Returns the minimum value of the distribution */ + IntType min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return _min; } + /** Returns the maximum value of the distribution */ + IntType max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return _max; } + + /** Returns the minimum value of the distribution */ + IntType a() const { return _min; } + /** Returns the maximum value of the distribution */ + IntType b() const { return _max; } + + /** Returns the parameters of the distribution. */ + param_type param() const { return param_type(_min, _max); } + /** Sets the parameters of the distribution. */ + void param(const param_type& parm) + { + _min = parm.a(); + _max = parm.b(); + } + + /** + * Effects: Subsequent uses of the distribution do not depend + * on values produced by any engine prior to invoking reset. + */ + void reset() { } + + /** Returns an integer uniformly distributed in the range [min, max]. */ + template<class Engine> + result_type operator()(Engine& eng) const + { return detail::generate_uniform_int(eng, _min, _max); } + + /** + * Returns an integer uniformly distributed in the range + * [param.a(), param.b()]. + */ + template<class Engine> + result_type operator()(Engine& eng, const param_type& parm) const + { return detail::generate_uniform_int(eng, parm.a(), parm.b()); } + + /** Writes the distribution to a @c std::ostream. */ + BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, uniform_int_distribution, ud) + { + os << ud.param(); + return os; + } + + /** Reads the distribution from a @c std::istream. */ + BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, uniform_int_distribution, ud) + { + param_type parm; + if(is >> parm) { + ud.param(parm); + } + return is; + } + + /** + * Returns true if the two distributions will produce identical sequences + * of values given equal generators. + */ + BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(uniform_int_distribution, lhs, rhs) + { return lhs._min == rhs._min && lhs._max == rhs._max; } + + /** + * Returns true if the two distributions may produce different sequences + * of values given equal generators. + */ + BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(uniform_int_distribution) + +private: + IntType _min; + IntType _max; +}; + +} // namespace random +} // namespace boost + +#endif // BOOST_RANDOM_UNIFORM_INT_HPP diff --git a/3rdParty/Boost/src/boost/random/variate_generator.hpp b/3rdParty/Boost/src/boost/random/variate_generator.hpp index 930d961..ac69800 100644 --- a/3rdParty/Boost/src/boost/random/variate_generator.hpp +++ b/3rdParty/Boost/src/boost/random/variate_generator.hpp @@ -1,34 +1,22 @@ /* boost random/variate_generator.hpp header file * * Copyright Jens Maurer 2002 + * Copyright Steven Watanabe 2011 * 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: variate_generator.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ + * $Id: variate_generator.hpp 71018 2011-04-05 21:27:52Z steven_watanabe $ * */ #ifndef BOOST_RANDOM_RANDOM_GENERATOR_HPP #define BOOST_RANDOM_RANDOM_GENERATOR_HPP -#include <boost/config.hpp> - -// implementation details -#include <boost/detail/workaround.hpp> -#include <boost/random/uniform_01.hpp> -#include <boost/random/detail/pass_through_engine.hpp> -#include <boost/random/detail/uniform_int_float.hpp> #include <boost/random/detail/ptr_helper.hpp> -// Borland C++ 5.6.0 has problems using its numeric_limits traits as -// template parameters -#if BOOST_WORKAROUND(__BORLANDC__, <= 0x564) -#include <boost/type_traits/is_integral.hpp> -#endif - #include <boost/random/detail/disable_warnings.hpp> namespace boost { @@ -36,54 +24,6 @@ namespace boost { /// \cond hide_private_members namespace random { -namespace detail { - -template<bool have_int, bool want_int> -struct engine_helper; - -// for consistency, always have two levels of decorations -template<> -struct engine_helper<true, true> -{ - template<class Engine, class DistInputType> - struct impl - { - typedef pass_through_engine<Engine> type; - }; -}; - -template<> -struct engine_helper<false, false> -{ - template<class Engine, class DistInputType> - struct impl - { - typedef uniform_01<Engine, DistInputType> type; - }; -}; - -template<> -struct engine_helper<true, false> -{ - template<class Engine, class DistInputType> - struct impl - { - typedef uniform_01<Engine, DistInputType> type; - }; -}; - -template<> -struct engine_helper<false, true> -{ - template<class Engine, class DistInputType> - struct impl - { - typedef uniform_int_float<Engine, unsigned long> type; - }; -}; - -} // namespace detail -} // namespace random ///\endcond @@ -93,9 +33,6 @@ struct engine_helper<false, true> * Boost.Random provides a vast choice of \generators as well * as \distributions. * - * Instantations of class template @c variate_generator model - * a \number_generator. - * * The argument for the template parameter Engine shall be of * the form U, U&, or U*, where U models a * \uniform_random_number_generator. Then, the member @@ -114,107 +51,72 @@ template<class Engine, class Distribution> class variate_generator { private: - typedef random::detail::pass_through_engine<Engine> decorated_engine; - + typedef boost::random::detail::ptr_helper<Engine> helper_type; public: - typedef typename decorated_engine::base_type engine_value_type; - typedef Engine engine_type; - typedef Distribution distribution_type; - typedef typename Distribution::result_type result_type; - - /** - * Constructs a @c variate_generator object with the associated - * \uniform_random_number_generator eng and the associated - * \random_distribution d. - * - * Throws: If and what the copy constructor of Engine or - * Distribution throws. - */ - variate_generator(Engine e, Distribution d) - : _eng(decorated_engine(e)), _dist(d) { } - - /** - * Returns: distribution()(e) - * - * Notes: The sequence of numbers produced by the - * \uniform_random_number_generator e, s<sub>e</sub>, is - * obtained from the sequence of numbers produced by the - * associated \uniform_random_number_generator eng, s<sub>eng</sub>, - * as follows: Consider the values of @c numeric_limits<T>::is_integer - * for @c T both @c Distribution::input_type and - * @c engine_value_type::result_type. If the values for both types are - * true, then se is identical to s<sub>eng</sub>. Otherwise, if the - * values for both types are false, then the numbers in s<sub>eng</sub> - * are divided by engine().max()-engine().min() to obtain the numbers - * in s<sub>e</sub>. Otherwise, if the value for - * @c engine_value_type::result_type is true and the value for - * @c Distribution::input_type is false, then the numbers in s<sub>eng</sub> - * are divided by engine().max()-engine().min()+1 to obtain the numbers in - * s<sub>e</sub>. Otherwise, the mapping from s<sub>eng</sub> to - * s<sub>e</sub> is implementation-defined. In all cases, an - * implicit conversion from @c engine_value_type::result_type to - * @c Distribution::input_type is performed. If such a conversion does - * not exist, the program is ill-formed. - */ - result_type operator()() { return _dist(_eng); } - /** - * Returns: distribution()(e, value). - * For the semantics of e, see the description of operator()(). - */ - template<class T> - result_type operator()(T value) { return _dist(_eng, value); } - - /** - * Returns: A reference to the associated uniform random number generator. - */ - engine_value_type& engine() { return _eng.base().base(); } - /** - * Returns: A reference to the associated uniform random number generator. - */ - const engine_value_type& engine() const { return _eng.base().base(); } - - /** - * Returns: A reference to the associated random distribution. - */ - distribution_type& distribution() { return _dist; } - /** - * Returns: A reference to the associated random distribution. - */ - const distribution_type& distribution() const { return _dist; } - - /** - * Precondition: distribution().min() is well-formed - * - * Returns: distribution().min() - */ - result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return (distribution().min)(); } - /** - * Precondition: distribution().max() is well-formed - * - * Returns: distribution().max() - */ - result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return (distribution().max)(); } + typedef typename helper_type::value_type engine_value_type; + typedef Engine engine_type; + typedef Distribution distribution_type; + typedef typename Distribution::result_type result_type; + + /** + * Constructs a @c variate_generator object with the associated + * \uniform_random_number_generator eng and the associated + * \random_distribution d. + * + * Throws: If and what the copy constructor of Engine or + * Distribution throws. + */ + variate_generator(Engine e, Distribution d) + : _eng(e), _dist(d) { } + + /** Returns: distribution()(engine()) */ + result_type operator()() { return _dist(engine()); } + /** + * Returns: distribution()(engine(), value). + */ + template<class T> + result_type operator()(const T& value) { return _dist(engine(), value); } + + /** + * Returns: A reference to the associated uniform random number generator. + */ + engine_value_type& engine() { return helper_type::ref(_eng); } + /** + * Returns: A reference to the associated uniform random number generator. + */ + const engine_value_type& engine() const { return helper_type::ref(_eng); } + + /** Returns: A reference to the associated \random_distribution. */ + distribution_type& distribution() { return _dist; } + /** + * Returns: A reference to the associated random distribution. + */ + const distribution_type& distribution() const { return _dist; } + + /** + * Precondition: distribution().min() is well-formed + * + * Returns: distribution().min() + */ + result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return (distribution().min)(); } + /** + * Precondition: distribution().max() is well-formed + * + * Returns: distribution().max() + */ + result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return (distribution().max)(); } private: -#if BOOST_WORKAROUND(__BORLANDC__, <= 0x564) - typedef typename random::detail::engine_helper< - ::boost::is_integral<typename decorated_engine::result_type>::value, - ::boost::is_integral<typename Distribution::input_type>::value - >::BOOST_NESTED_TEMPLATE impl<decorated_engine, typename Distribution::input_type>::type internal_engine_type; -#else - enum { - have_int = std::numeric_limits<typename decorated_engine::result_type>::is_integer, - want_int = std::numeric_limits<typename Distribution::input_type>::is_integer - }; - typedef typename random::detail::engine_helper<have_int, want_int>::BOOST_NESTED_TEMPLATE impl<decorated_engine, typename Distribution::input_type>::type internal_engine_type; -#endif - - internal_engine_type _eng; - distribution_type _dist; + Engine _eng; + distribution_type _dist; }; +} // namespace random + +using random::variate_generator; + } // namespace boost -#include <boost/random/detail/disable_warnings.hpp> +#include <boost/random/detail/enable_warnings.hpp> #endif // BOOST_RANDOM_RANDOM_GENERATOR_HPP |