diff options
author | Stefan Radomski <github@mintwerk.de> | 2016-05-12 13:12:33 (GMT) |
---|---|---|
committer | Stefan Radomski <github@mintwerk.de> | 2016-05-12 13:12:33 (GMT) |
commit | b62e7979600feee23dc7cdb61042a8fc7673122b (patch) | |
tree | f7351372f37979dd2d048e0b68a16a4cd3b2aadb /contrib/src/boost/random/detail | |
parent | 1b11b310be61e51b3ac5ebb83f7c8a33aef3d6e8 (diff) | |
download | uscxml-b62e7979600feee23dc7cdb61042a8fc7673122b.zip uscxml-b62e7979600feee23dc7cdb61042a8fc7673122b.tar.gz uscxml-b62e7979600feee23dc7cdb61042a8fc7673122b.tar.bz2 |
Major Refactoring v2.0
Diffstat (limited to 'contrib/src/boost/random/detail')
-rw-r--r-- | contrib/src/boost/random/detail/config.hpp | 18 | ||||
-rw-r--r-- | contrib/src/boost/random/detail/const_mod.hpp | 216 | ||||
-rw-r--r-- | contrib/src/boost/random/detail/disable_warnings.hpp | 29 | ||||
-rw-r--r-- | contrib/src/boost/random/detail/enable_warnings.hpp | 22 | ||||
-rw-r--r-- | contrib/src/boost/random/detail/generator_bits.hpp | 36 | ||||
-rw-r--r-- | contrib/src/boost/random/detail/generator_seed_seq.hpp | 40 | ||||
-rw-r--r-- | contrib/src/boost/random/detail/integer_log2.hpp | 84 | ||||
-rw-r--r-- | contrib/src/boost/random/detail/large_arithmetic.hpp | 122 | ||||
-rw-r--r-- | contrib/src/boost/random/detail/operators.hpp | 84 | ||||
-rw-r--r-- | contrib/src/boost/random/detail/polynomial.hpp | 384 | ||||
-rw-r--r-- | contrib/src/boost/random/detail/ptr_helper.hpp | 67 | ||||
-rw-r--r-- | contrib/src/boost/random/detail/seed.hpp | 115 | ||||
-rw-r--r-- | contrib/src/boost/random/detail/seed_impl.hpp | 398 | ||||
-rw-r--r-- | contrib/src/boost/random/detail/signed_unsigned_tools.hpp | 89 | ||||
-rw-r--r-- | contrib/src/boost/random/detail/uniform_int_float.hpp | 76 |
15 files changed, 1780 insertions, 0 deletions
diff --git a/contrib/src/boost/random/detail/config.hpp b/contrib/src/boost/random/detail/config.hpp new file mode 100644 index 0000000..724ab19 --- /dev/null +++ b/contrib/src/boost/random/detail/config.hpp @@ -0,0 +1,18 @@ +/* boost random/detail/config.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$ + */ + +#include <boost/config.hpp> + +#if (defined(BOOST_NO_OPERATORS_IN_NAMESPACE) || defined(BOOST_NO_MEMBER_TEMPLATE_FRIENDS)) \ + && !defined(BOOST_MSVC) + #define BOOST_RANDOM_NO_STREAM_OPERATORS +#endif diff --git a/contrib/src/boost/random/detail/const_mod.hpp b/contrib/src/boost/random/detail/const_mod.hpp new file mode 100644 index 0000000..e0a43ab --- /dev/null +++ b/contrib/src/boost/random/detail/const_mod.hpp @@ -0,0 +1,216 @@ +/* boost random/detail/const_mod.hpp header file + * + * Copyright Jens Maurer 2000-2001 + * 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$ + * + * Revision history + * 2001-02-18 moved to individual header files + */ + +#ifndef BOOST_RANDOM_CONST_MOD_HPP +#define BOOST_RANDOM_CONST_MOD_HPP + +#include <boost/assert.hpp> +#include <boost/static_assert.hpp> +#include <boost/integer_traits.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 { + +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 suppress_warnings = (m == 0); + BOOST_ASSERT(suppress_warnings == 0); + return x % (m + suppress_warnings); + } + } + + static IntType add(IntType x, IntType c) + { + 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(x < m - c) + return x + c; + else + return x - (m - c); + } + + static IntType mult(IntType a, IntType x) + { + 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 + return mult_general(a, x); + } + + static IntType mult_add(IntType a, IntType x, IntType c) + { + 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 suppress_warnings = (m == 0); + BOOST_ASSERT(suppress_warnings == 0); + return (a*x+c) % (m + suppress_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 : (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 mult_small(IntType a, IntType x) + { + IntType suppress_warnings = (m == 0); + BOOST_ASSERT(suppress_warnings == 0); + return a*x % (m + suppress_warnings); + } + + static IntType mult_schrage(IntType a, IntType value) + { + const IntType q = m / a; + const IntType r = m % a; + + BOOST_ASSERT(r < q); // check that overflow cannot happen + + 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)); + } + } + + 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_ASSERT(c > 0); + IntType l1 = 0; + IntType l2 = 1; + IntType n = c; + IntType p = m; + for(;;) { + IntType q = p / n; + l1 += q * l2; + p -= q * n; + if(p == 0) + return l2; + IntType q2 = n / p; + l2 += q2 * l1; + n -= q2 * p; + if(n == 0) + return m - l1; + } + } + + // 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_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(;;) { + if(p == 0) + return l2; + IntType q2 = n / p; + l2 += q2 * l1; + n -= q2 * p; + if(n == 0) + return m - l1; + q = p / n; + l1 += q * l2; + p -= q * n; + } + } +}; + +} // namespace random +} // namespace boost + +#include <boost/random/detail/enable_warnings.hpp> + +#endif // BOOST_RANDOM_CONST_MOD_HPP diff --git a/contrib/src/boost/random/detail/disable_warnings.hpp b/contrib/src/boost/random/detail/disable_warnings.hpp new file mode 100644 index 0000000..4582dcb --- /dev/null +++ b/contrib/src/boost/random/detail/disable_warnings.hpp @@ -0,0 +1,29 @@ +/* boost random/detail/disable_warnings.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$ + * + */ + +// No #include guard. This header is intended to be included multiple times. + +#include <boost/config.hpp> + +#ifdef BOOST_MSVC +#pragma warning(push) +#pragma warning(disable:4512) +#pragma warning(disable:4127) +#pragma warning(disable:4724) +#pragma warning(disable:4800) // 'int' : forcing value to bool 'true' or 'false' (performance warning) +#endif + +#if defined(BOOST_GCC) && BOOST_GCC >= 40600 +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wlogical-op" +#endif diff --git a/contrib/src/boost/random/detail/enable_warnings.hpp b/contrib/src/boost/random/detail/enable_warnings.hpp new file mode 100644 index 0000000..24f3bb3 --- /dev/null +++ b/contrib/src/boost/random/detail/enable_warnings.hpp @@ -0,0 +1,22 @@ +/* boost random/detail/enable_warnings.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$ + * + */ + +// No #include guard. This header is intended to be included multiple times. + +#ifdef BOOST_MSVC +#pragma warning(pop) +#endif + +#if defined(BOOST_GCC) && BOOST_GCC >= 40600 +#pragma GCC diagnostic pop +#endif diff --git a/contrib/src/boost/random/detail/generator_bits.hpp b/contrib/src/boost/random/detail/generator_bits.hpp new file mode 100644 index 0000000..0527614 --- /dev/null +++ b/contrib/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$ + * + */ + +#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/contrib/src/boost/random/detail/generator_seed_seq.hpp b/contrib/src/boost/random/detail/generator_seed_seq.hpp new file mode 100644 index 0000000..7e13483 --- /dev/null +++ b/contrib/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$ + * + */ + +#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/contrib/src/boost/random/detail/integer_log2.hpp b/contrib/src/boost/random/detail/integer_log2.hpp new file mode 100644 index 0000000..248243a --- /dev/null +++ b/contrib/src/boost/random/detail/integer_log2.hpp @@ -0,0 +1,84 @@ +/* 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$ + * + */ + +#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 { + +#if !defined(BOOST_NO_CXX11_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 inline __attribute__((__const__)) __attribute__((__always_inline__)) +#else +#define BOOST_RANDOM_DETAIL_CONSTEXPR inline +#endif + +template<int Shift> +struct integer_log2_impl +{ +#if defined(BOOST_NO_CXX11_CONSTEXPR) + 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); + } +#else + template<class T> + BOOST_RANDOM_DETAIL_CONSTEXPR static int apply2(T t, int accum, int update) + { + return integer_log2_impl<Shift / 2>::apply(t >> update, accum + update); + } + + template<class T> + BOOST_RANDOM_DETAIL_CONSTEXPR static int apply(T t, int accum) + { + return apply2(t, accum, ((t >> Shift) != 0) * Shift); + } +#endif +}; + +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/contrib/src/boost/random/detail/large_arithmetic.hpp b/contrib/src/boost/random/detail/large_arithmetic.hpp new file mode 100644 index 0000000..66f6b4e --- /dev/null +++ b/contrib/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$ + */ + +#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) +{ + const int bits = + ::std::numeric_limits< ::boost::uintmax_t>::digits / 2; + 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] = static_cast<digit_t>((rem >> bits) & mask); + product[i-2] = static_cast<digit_t>(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/contrib/src/boost/random/detail/operators.hpp b/contrib/src/boost/random/detail/operators.hpp new file mode 100644 index 0000000..597343c --- /dev/null +++ b/contrib/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$ + */ + +#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/contrib/src/boost/random/detail/polynomial.hpp b/contrib/src/boost/random/detail/polynomial.hpp new file mode 100644 index 0000000..a8c4b26 --- /dev/null +++ b/contrib/src/boost/random/detail/polynomial.hpp @@ -0,0 +1,384 @@ +/* boost random/detail/polynomial.hpp header file + * + * Copyright Steven Watanabe 2014 + * Distributed under the Boost Software License, Version 1.0. (See + * accompanying file LICENSE_1_0.txt or copy at + * http://www.boost.org/LICENSE_1_0.txt) + * + * See http://www.boost.org for most recent version including documentation. + * + * $Id$ + */ + +#ifndef BOOST_RANDOM_DETAIL_POLYNOMIAL_HPP +#define BOOST_RANDOM_DETAIL_POLYNOMIAL_HPP + +#include <cstddef> +#include <limits> +#include <vector> +#include <algorithm> +#include <boost/assert.hpp> +#include <boost/cstdint.hpp> + +namespace boost { +namespace random { +namespace detail { + +class polynomial_ops { +public: + typedef unsigned long digit_t; + + static void add(std::size_t size, const digit_t * lhs, + const digit_t * rhs, digit_t * output) + { + for(std::size_t i = 0; i < size; ++i) { + output[i] = lhs[i] ^ rhs[i]; + } + } + + static void add_shifted_inplace(std::size_t size, const digit_t * lhs, + digit_t * output, std::size_t shift) + { + if(shift == 0) { + add(size, lhs, output, output); + return; + } + std::size_t bits = std::numeric_limits<digit_t>::digits; + digit_t prev = 0; + for(std::size_t i = 0; i < size; ++i) { + digit_t tmp = lhs[i]; + output[i] ^= (tmp << shift) | (prev >> (bits-shift)); + prev = tmp; + } + output[size] ^= (prev >> (bits-shift)); + } + + static void multiply_simple(std::size_t size, const digit_t * lhs, + const digit_t * rhs, digit_t * output) + { + std::size_t bits = std::numeric_limits<digit_t>::digits; + for(std::size_t i = 0; i < 2*size; ++i) { + output[i] = 0; + } + for(std::size_t i = 0; i < size; ++i) { + for(std::size_t j = 0; j < bits; ++j) { + if((lhs[i] & (digit_t(1) << j)) != 0) { + add_shifted_inplace(size, rhs, output + i, j); + } + } + } + } + + // memory requirements: (size - cutoff) * 4 + next_smaller + static void multiply_karatsuba(std::size_t size, + const digit_t * lhs, const digit_t * rhs, + digit_t * output) + { + if(size < 64) { + multiply_simple(size, lhs, rhs, output); + return; + } + // split in half + std::size_t cutoff = size/2; + multiply_karatsuba(cutoff, lhs, rhs, output); + multiply_karatsuba(size - cutoff, lhs + cutoff, rhs + cutoff, + output + cutoff*2); + std::vector<digit_t> local1(size - cutoff); + std::vector<digit_t> local2(size - cutoff); + // combine the digits for the inner multiply + add(cutoff, lhs, lhs + cutoff, &local1[0]); + if(size & 1) local1[cutoff] = lhs[size - 1]; + add(cutoff, rhs + cutoff, rhs, &local2[0]); + if(size & 1) local2[cutoff] = rhs[size - 1]; + std::vector<digit_t> local3((size - cutoff) * 2); + multiply_karatsuba(size - cutoff, &local1[0], &local2[0], &local3[0]); + add(cutoff * 2, output, &local3[0], &local3[0]); + add((size - cutoff) * 2, output + cutoff*2, &local3[0], &local3[0]); + // Finally, add the inner result + add((size - cutoff) * 2, output + cutoff, &local3[0], output + cutoff); + } + + static void multiply_add_karatsuba(std::size_t size, + const digit_t * lhs, const digit_t * rhs, + digit_t * output) + { + std::vector<digit_t> buf(size * 2); + multiply_karatsuba(size, lhs, rhs, &buf[0]); + add(size * 2, &buf[0], output, output); + } + + static void multiply(const digit_t * lhs, std::size_t lhs_size, + const digit_t * rhs, std::size_t rhs_size, + digit_t * output) + { + std::fill_n(output, lhs_size + rhs_size, digit_t(0)); + multiply_add(lhs, lhs_size, rhs, rhs_size, output); + } + + static void multiply_add(const digit_t * lhs, std::size_t lhs_size, + const digit_t * rhs, std::size_t rhs_size, + digit_t * output) + { + // split into pieces that can be passed to + // karatsuba multiply. + while(lhs_size != 0) { + if(lhs_size < rhs_size) { + std::swap(lhs, rhs); + std::swap(lhs_size, rhs_size); + } + + multiply_add_karatsuba(rhs_size, lhs, rhs, output); + + lhs += rhs_size; + lhs_size -= rhs_size; + output += rhs_size; + } + } + + static void copy_bits(const digit_t * x, std::size_t low, std::size_t high, + digit_t * out) + { + const std::size_t bits = std::numeric_limits<digit_t>::digits; + std::size_t offset = low/bits; + x += offset; + low -= offset*bits; + high -= offset*bits; + std::size_t n = (high-low)/bits; + if(low == 0) { + for(std::size_t i = 0; i < n; ++i) { + out[i] = x[i]; + } + } else { + for(std::size_t i = 0; i < n; ++i) { + out[i] = (x[i] >> low) | (x[i+1] << (bits-low)); + } + } + if((high-low)%bits) { + digit_t low_mask = (digit_t(1) << ((high-low)%bits)) - 1; + digit_t result = (x[n] >> low); + if(low != 0 && (n+1)*bits < high) { + result |= (x[n+1] << (bits-low)); + } + out[n] = (result & low_mask); + } + } + + static void shift_left(digit_t * val, std::size_t size, std::size_t shift) + { + const std::size_t bits = std::numeric_limits<digit_t>::digits; + BOOST_ASSERT(shift > 0); + BOOST_ASSERT(shift < bits); + digit_t prev = 0; + for(std::size_t i = 0; i < size; ++i) { + digit_t tmp = val[i]; + val[i] = (prev >> (bits - shift)) | (val[i] << shift); + prev = tmp; + } + } + + static digit_t sqr(digit_t val) { + const std::size_t bits = std::numeric_limits<digit_t>::digits; + digit_t mask = (digit_t(1) << bits/2) - 1; + for(std::size_t i = bits; i > 1; i /= 2) { + val = ((val & ~mask) << i/2) | (val & mask); + mask = mask & (mask >> i/4); + mask = mask | (mask << i/2); + } + return val; + } + + static void sqr(digit_t * val, std::size_t size) + { + const std::size_t bits = std::numeric_limits<digit_t>::digits; + digit_t mask = (digit_t(1) << bits/2) - 1; + for(std::size_t i = 0; i < size; ++i) { + digit_t x = val[size - i - 1]; + val[(size - i - 1) * 2] = sqr(x & mask); + val[(size - i - 1) * 2 + 1] = sqr(x >> bits/2); + } + } + + // optimized for the case when the modulus has few bits set. + struct sparse_mod { + sparse_mod(const digit_t * divisor, std::size_t divisor_bits) + { + const std::size_t bits = std::numeric_limits<digit_t>::digits; + _remainder_bits = divisor_bits - 1; + for(std::size_t i = 0; i < divisor_bits; ++i) { + if(divisor[i/bits] & (digit_t(1) << i%bits)) { + _bit_indices.push_back(i); + } + } + BOOST_ASSERT(_bit_indices.back() == divisor_bits - 1); + _bit_indices.pop_back(); + if(_bit_indices.empty()) { + _block_bits = divisor_bits; + _lower_bits = 0; + } else { + _block_bits = divisor_bits - _bit_indices.back() - 1; + _lower_bits = _bit_indices.back() + 1; + } + + _partial_quotient.resize((_block_bits + bits - 1)/bits); + } + void operator()(digit_t * dividend, std::size_t dividend_bits) + { + const std::size_t bits = std::numeric_limits<digit_t>::digits; + while(dividend_bits > _remainder_bits) { + std::size_t block_start = (std::max)(dividend_bits - _block_bits, _remainder_bits); + std::size_t block_size = (dividend_bits - block_start + bits - 1) / bits; + copy_bits(dividend, block_start, dividend_bits, &_partial_quotient[0]); + for(std::size_t i = 0; i < _bit_indices.size(); ++i) { + std::size_t pos = _bit_indices[i] + block_start - _remainder_bits; + add_shifted_inplace(block_size, &_partial_quotient[0], dividend + pos/bits, pos%bits); + } + add_shifted_inplace(block_size, &_partial_quotient[0], dividend + block_start/bits, block_start%bits); + dividend_bits = block_start; + } + } + std::vector<digit_t> _partial_quotient; + std::size_t _remainder_bits; + std::size_t _block_bits; + std::size_t _lower_bits; + std::vector<std::size_t> _bit_indices; + }; + + // base should have the same number of bits as mod + // base, and mod should both be able to hold a power + // of 2 >= mod_bits. out needs to be twice as large. + static void mod_pow_x(boost::uintmax_t exponent, const digit_t * mod, std::size_t mod_bits, digit_t * out) + { + const std::size_t bits = std::numeric_limits<digit_t>::digits; + const std::size_t n = (mod_bits + bits - 1) / bits; + const std::size_t highbit = mod_bits - 1; + if(exponent == 0) { + out[0] = 1; + std::fill_n(out + 1, n - 1, digit_t(0)); + return; + } + boost::uintmax_t i = std::numeric_limits<boost::uintmax_t>::digits - 1; + while(((boost::uintmax_t(1) << i) & exponent) == 0) { + --i; + } + out[0] = 2; + std::fill_n(out + 1, n - 1, digit_t(0)); + sparse_mod m(mod, mod_bits); + while(i--) { + sqr(out, n); + m(out, 2 * mod_bits - 1); + if((boost::uintmax_t(1) << i) & exponent) { + shift_left(out, n, 1); + if(out[highbit / bits] & (digit_t(1) << highbit%bits)) + add(n, out, mod, out); + } + } + } +}; + +class polynomial +{ + typedef polynomial_ops::digit_t digit_t; +public: + polynomial() : _size(0) {} + class reference { + public: + reference(digit_t &value, int idx) + : _value(value), _idx(idx) {} + operator bool() const { return (_value & (digit_t(1) << _idx)) != 0; } + reference& operator=(bool b) + { + if(b) { + _value |= (digit_t(1) << _idx); + } else { + _value &= ~(digit_t(1) << _idx); + } + return *this; + } + reference &operator^=(bool b) + { + _value ^= (digit_t(b) << _idx); + return *this; + } + + reference &operator=(const reference &other) + { + return *this = static_cast<bool>(other); + } + private: + digit_t &_value; + int _idx; + }; + reference operator[](std::size_t i) + { + static const std::size_t bits = std::numeric_limits<digit_t>::digits; + ensure_bit(i); + return reference(_storage[i/bits], i%bits); + } + bool operator[](std::size_t i) const + { + static const std::size_t bits = std::numeric_limits<digit_t>::digits; + if(i < size()) + return (_storage[i/bits] & (digit_t(1) << (i%bits))) != 0; + else + return false; + } + std::size_t size() const + { + return _size; + } + void resize(std::size_t n) + { + static const std::size_t bits = std::numeric_limits<digit_t>::digits; + _storage.resize((n + bits - 1)/bits); + // clear the high order bits in case we're shrinking. + if(n%bits) { + _storage.back() &= ((digit_t(1) << (n%bits)) - 1); + } + _size = n; + } + friend polynomial operator*(const polynomial &lhs, const polynomial &rhs); + friend polynomial mod_pow_x(boost::uintmax_t exponent, polynomial mod); +private: + std::vector<polynomial_ops::digit_t> _storage; + std::size_t _size; + void ensure_bit(std::size_t i) + { + if(i >= size()) { + resize(i + 1); + } + } + void normalize() + { + while(size() && (*this)[size() - 1] == 0) + resize(size() - 1); + } +}; + +inline polynomial operator*(const polynomial &lhs, const polynomial &rhs) +{ + polynomial result; + result._storage.resize(lhs._storage.size() + rhs._storage.size()); + polynomial_ops::multiply(&lhs._storage[0], lhs._storage.size(), + &rhs._storage[0], rhs._storage.size(), + &result._storage[0]); + result._size = lhs._size + rhs._size; + return result; +} + +inline polynomial mod_pow_x(boost::uintmax_t exponent, polynomial mod) +{ + polynomial result; + mod.normalize(); + std::size_t mod_size = mod.size(); + result._storage.resize(mod._storage.size() * 2); + result._size = mod.size() * 2; + polynomial_ops::mod_pow_x(exponent, &mod._storage[0], mod_size, &result._storage[0]); + result.resize(mod.size() - 1); + return result; +} + +} +} +} + +#endif // BOOST_RANDOM_DETAIL_POLYNOMIAL_HPP diff --git a/contrib/src/boost/random/detail/ptr_helper.hpp b/contrib/src/boost/random/detail/ptr_helper.hpp new file mode 100644 index 0000000..f1b983d --- /dev/null +++ b/contrib/src/boost/random/detail/ptr_helper.hpp @@ -0,0 +1,67 @@ +/* boost random/detail/ptr_helper.hpp header file + * + * Copyright Jens Maurer 2002 + * Distributed under the Boost Software License, Version 1.0. (See + * accompanying file LICENSE_1_0.txt or copy at + * http://www.boost.org/LICENSE_1_0.txt) + * + * See http://www.boost.org for most recent version including documentation. + * + * $Id$ + * + */ + +#ifndef BOOST_RANDOM_DETAIL_PTR_HELPER_HPP +#define BOOST_RANDOM_DETAIL_PTR_HELPER_HPP + +#include <boost/config.hpp> + + +namespace boost { +namespace random { +namespace detail { + +// type_traits could help here, but I don't want to depend on type_traits. +template<class T> +struct ptr_helper +{ + typedef T value_type; + typedef T& reference_type; + typedef const T& rvalue_type; + static reference_type ref(T& r) { return r; } + static const T& ref(const T& r) { return r; } +}; + +template<class T> +struct ptr_helper<T&> +{ + typedef T value_type; + typedef T& reference_type; + typedef T& rvalue_type; + static reference_type ref(T& r) { return r; } + static const T& ref(const T& r) { return r; } +}; + +template<class T> +struct ptr_helper<T*> +{ + typedef T value_type; + typedef T& reference_type; + typedef T* rvalue_type; + static reference_type ref(T * p) { return *p; } + static const T& ref(const T * p) { return *p; } +}; + +} // namespace detail +} // namespace random +} // namespace boost + +// +// BOOST_RANDOM_PTR_HELPER_SPEC -- +// +// Helper macro for broken compilers defines specializations of +// ptr_helper. +// +# define BOOST_RANDOM_PTR_HELPER_SPEC(T) + +#endif // BOOST_RANDOM_DETAIL_PTR_HELPER_HPP diff --git a/contrib/src/boost/random/detail/seed.hpp b/contrib/src/boost/random/detail/seed.hpp new file mode 100644 index 0000000..557482a --- /dev/null +++ b/contrib/src/boost/random/detail/seed.hpp @@ -0,0 +1,115 @@ +/* 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$ + */ + +#ifndef BOOST_RANDOM_DETAIL_SEED_HPP +#define BOOST_RANDOM_DETAIL_SEED_HPP + +#include <boost/config.hpp> + +// Sun seems to have trouble with the use of SFINAE for the +// templated constructor. So does Borland. +#if !defined(BOOST_NO_SFINAE) && !defined(__SUNPRO_CC) && !defined(__BORLANDC__) + +#include <boost/utility/enable_if.hpp> +#include <boost/type_traits/is_arithmetic.hpp> +#include <boost/mpl/bool.hpp> + +namespace boost { +namespace random { +namespace detail { + +template<class T> +struct disable_seed : boost::disable_if<boost::is_arithmetic<T> > {}; + +template<class Engine, class T> +struct disable_constructor : disable_seed<T> {}; + +template<class Engine> +struct disable_constructor<Engine, Engine> {}; + +#define BOOST_RANDOM_DETAIL_GENERATOR_CONSTRUCTOR(Self, Generator, gen) \ + template<class Generator> \ + explicit Self(Generator& gen, typename ::boost::random::detail::disable_constructor<Self, Generator>::type* = 0) + +#define BOOST_RANDOM_DETAIL_GENERATOR_SEED(Self, Generator, gen) \ + 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) +} +} +} + +#else + +#include <boost/type_traits/is_arithmetic.hpp> +#include <boost/mpl/bool.hpp> + +#define BOOST_RANDOM_DETAIL_GENERATOR_CONSTRUCTOR(Self, Generator, gen) \ + Self(Self& other) { *this = other; } \ + Self(const Self& other) { *this = other; } \ + template<class Generator> \ + explicit Self(Generator& gen) { \ + boost_random_constructor_impl(gen, ::boost::is_arithmetic<Generator>());\ + } \ + template<class Generator> \ + void boost_random_constructor_impl(Generator& gen, ::boost::mpl::false_) + +#define BOOST_RANDOM_DETAIL_GENERATOR_SEED(Self, Generator, gen) \ + template<class Generator> \ + void seed(Generator& gen) { \ + boost_random_seed_impl(gen, ::boost::is_arithmetic<Generator>());\ + }\ + 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_) + +#define BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(Self, T, x) \ + void seed(const T& x) { boost_random_seed_impl(x, ::boost::mpl::true_()); }\ + void boost_random_seed_impl(const T& x, ::boost::mpl::true_) + +#endif + +#endif diff --git a/contrib/src/boost/random/detail/seed_impl.hpp b/contrib/src/boost/random/detail/seed_impl.hpp new file mode 100644 index 0000000..918a294 --- /dev/null +++ b/contrib/src/boost/random/detail/seed_impl.hpp @@ -0,0 +1,398 @@ +/* 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$ + */ + +#ifndef BOOST_RANDOM_DETAIL_SEED_IMPL_HPP +#define BOOST_RANDOM_DETAIL_SEED_IMPL_HPP + +#include <stdexcept> +#include <boost/cstdint.hpp> +#include <boost/throw_exception.hpp> +#include <boost/config/no_tr1/cmath.hpp> +#include <boost/integer/integer_mask.hpp> +#include <boost/integer/static_log2.hpp> +#include <boost/random/traits.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(T(arg * arg), n / 2, + n%2 == 0? result : T(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::random::traits::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 suppress_warning = (bits >= 32); + BOOST_ASSERT(suppress_warning == 1); + val >>= (32 * suppress_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 suppress_warning = (bits >= 32); + BOOST_ASSERT(suppress_warning == 1); + val >>= (32 * suppress_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::random::traits::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) { + boost::throw_exception(::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_) +{ + BOOST_STATIC_ASSERT_MSG(boost::is_integral<IntType>::value, "Sorry but this routine has not been ported to non built-in integers as it relies on a reinterpret_cast."); + 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::random::traits::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) { + boost::throw_exception(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_) +{ + BOOST_STATIC_ASSERT_MSG(boost::is_integral<IntType>::value, "Sorry but this routine has not been ported to non built-in integers as it relies on a reinterpret_cast."); + 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::random::traits::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) boost::throw_exception(std::invalid_argument("Not enough elements in call to seed.")); + val += *first * mult; + mult *= two32; + } + if(mask != 0) { + if(first == last) boost::throw_exception(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/contrib/src/boost/random/detail/signed_unsigned_tools.hpp b/contrib/src/boost/random/detail/signed_unsigned_tools.hpp new file mode 100644 index 0000000..1979908 --- /dev/null +++ b/contrib/src/boost/random/detail/signed_unsigned_tools.hpp @@ -0,0 +1,89 @@ +/* boost random/detail/signed_unsigned_tools.hpp header file + * + * Copyright Jens Maurer 2006 + * 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. + */ + +#ifndef BOOST_RANDOM_DETAIL_SIGNED_UNSIGNED_TOOLS +#define BOOST_RANDOM_DETAIL_SIGNED_UNSIGNED_TOOLS + +#include <boost/limits.hpp> +#include <boost/config.hpp> +#include <boost/random/traits.hpp> + +namespace boost { +namespace random { +namespace detail { + + +/* + * Compute x - y, we know that x >= y, return an unsigned value. + */ + +template<class T, bool sgn = std::numeric_limits<T>::is_signed && std::numeric_limits<T>::is_bounded> +struct subtract { }; + +template<class T> +struct subtract<T, /* signed */ false> +{ + typedef T result_type; + result_type operator()(T x, T y) { return x - y; } +}; + +template<class T> +struct subtract<T, /* signed */ true> +{ + typedef typename boost::random::traits::make_unsigned_or_unbounded<T>::type result_type; + result_type operator()(T x, T y) + { + if (y >= 0) // because x >= y, it follows that x >= 0, too + return result_type(x) - result_type(y); + if (x >= 0) // y < 0 + // avoid the nasty two's complement case for y == min() + return result_type(x) + result_type(-(y+1)) + 1; + // both x and y are negative: no signed overflow + return result_type(x - y); + } +}; + +/* + * Compute x + y, x is unsigned, result fits in type of "y". + */ + +template<class T1, class T2, bool sgn = (std::numeric_limits<T2>::is_signed && (std::numeric_limits<T1>::digits >= std::numeric_limits<T2>::digits))> +struct add { }; + +template<class T1, class T2> +struct add<T1, T2, /* signed or else T2 has more digits than T1 so the cast always works - needed when T2 is a multiprecision type and T1 is a native integer */ false> +{ + typedef T2 result_type; + result_type operator()(T1 x, T2 y) { return T2(x) + y; } +}; + +template<class T1, class T2> +struct add<T1, T2, /* signed */ true> +{ + typedef T2 result_type; + result_type operator()(T1 x, T2 y) + { + if (y >= 0) + return T2(x) + y; + // y < 0 + if (x > T1(-(y+1))) // result >= 0 after subtraction + // avoid the nasty two's complement edge case for y == min() + return T2(x - T1(-(y+1)) - 1); + // abs(x) < abs(y), thus T2 able to represent x + return T2(x) + y; + } +}; + +} // namespace detail +} // namespace random +} // namespace boost + +#endif // BOOST_RANDOM_DETAIL_SIGNED_UNSIGNED_TOOLS + diff --git a/contrib/src/boost/random/detail/uniform_int_float.hpp b/contrib/src/boost/random/detail/uniform_int_float.hpp new file mode 100644 index 0000000..393c455 --- /dev/null +++ b/contrib/src/boost/random/detail/uniform_int_float.hpp @@ -0,0 +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$ + * + */ + +#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/detail/generator_bits.hpp> + +#include <boost/random/detail/disable_warnings.hpp> + +namespace boost { +namespace random { +namespace detail { + +template<class URNG> +class uniform_int_float +{ +public: + typedef URNG base_type; + typedef typename base_type::result_type base_result; + + 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; + + uniform_int_float(base_type& rng) + : _rng(rng) {} + + 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; } + + result_type operator()() + { + base_result range = static_cast<base_result>((max)())+1; + return static_cast<result_type>(_rng() * range); + } + +private: + base_type& _rng; +}; + +} // namespace detail +} // namespace random +} // namespace boost + +#include <boost/random/detail/enable_warnings.hpp> + +#endif // BOOST_RANDOM_DETAIL_UNIFORM_INT_FLOAT_HPP |