diff --git a/sprout/detail/integer/static_log2.hpp b/sprout/detail/integer/static_log2.hpp new file mode 100644 index 00000000..67923ab7 --- /dev/null +++ b/sprout/detail/integer/static_log2.hpp @@ -0,0 +1,52 @@ +#ifndef SPROUT_DETAIL_INTEGER_STATIC_LOG2_HPP +#define SPROUT_DETAIL_INTEGER_STATIC_LOG2_HPP + +#include +#include + +namespace sprout { + namespace detail { + typedef std::uintmax_t static_log2_argument_type; + typedef int static_log2_result_type; + + namespace static_log2_impl { + typedef sprout::detail::static_log2_argument_type argument_type; + typedef sprout::detail::static_log2_result_type result_type; + + template + struct choose_initial_n { + SPROUT_STATIC_CONSTEXPR bool c = (sprout::detail::static_log2_impl::argument_type(1) << n << n) != 0; + SPROUT_STATIC_CONSTEXPR sprout::detail::static_log2_impl::result_type value = !c * n + choose_initial_n<2 * c * n>::value; + }; + template<> + struct choose_initial_n<0> { + SPROUT_STATIC_CONSTEXPR sprout::detail::static_log2_impl::result_type value = 0; + }; + + SPROUT_STATIC_CONSTEXPR sprout::detail::static_log2_impl::result_type n_zero = 16; + SPROUT_STATIC_CONSTEXPR sprout::detail::static_log2_impl::result_type initial_n = sprout::detail::static_log2_impl::choose_initial_n::value; + + template< + sprout::detail::static_log2_impl::argument_type x, + sprout::detail::static_log2_impl::result_type n = sprout::detail::static_log2_impl::initial_n + > + struct static_log2_impl { + SPROUT_STATIC_CONSTEXPR bool c = (x >> n) > 0; + SPROUT_STATIC_CONSTEXPR sprout::detail::static_log2_impl::result_type value = c * n + (static_log2_impl<(x >> c * n), n / 2>::value); + }; + template<> + struct static_log2_impl<1, 0> { + SPROUT_STATIC_CONSTEXPR sprout::detail::static_log2_impl::result_type value = 0; + }; + } // namespace static_log2_impl + + template + struct static_log2 { + SPROUT_STATIC_CONSTEXPR sprout::detail::static_log2_result_type value = sprout::detail::static_log2_impl::static_log2_impl::value; + }; + template<> + struct static_log2<0> {}; + } // namespace detail +} // namespace sprout + +#endif // #ifndef SPROUT_DETAIL_INTEGER_STATIC_LOG2_HPP diff --git a/sprout/random.hpp b/sprout/random.hpp index 8baba2e5..f24457d6 100644 --- a/sprout/random.hpp +++ b/sprout/random.hpp @@ -3,7 +3,9 @@ #include #include +#include #include +#include #include #include #include diff --git a/sprout/random/additive_combine.hpp b/sprout/random/additive_combine.hpp new file mode 100644 index 00000000..8db8a8f5 --- /dev/null +++ b/sprout/random/additive_combine.hpp @@ -0,0 +1,108 @@ +#ifndef SPROUT_RANDOM_ADDITIVE_COMBINE_HPP +#define SPROUT_RANDOM_ADDITIVE_COMBINE_HPP + +#include +#include +#include +#include +#include +#include + +namespace sprout { + namespace random { + // + // additive_combine_engine + // + template + class additive_combine_engine { + public: + typedef MLCG1 first_base; + typedef MLCG2 second_base; + typedef typename first_base::result_type result_type; + private: + struct private_constructor_tag {}; + private: + first_base mlcg1_; + second_base mlcg2_; + private: + SPROUT_CONSTEXPR additive_combine_engine( + first_base const& mlcg1, + second_base const& mlcg2, + private_constructor_tag + ) + : mlcg1_(mlcg1) + , mlcg2_(mlcg2) + {} + template + SPROUT_CONSTEXPR sprout::random::random_result generate(Random1 const& rnd1, Random2 const& rnd2) const { + return sprout::random::random_result( + rnd2.result() < rnd1.result() + ? rnd1.result() - rnd2.result() + : rnd1.result() - rnd2.result() + first_base::modulus - 1 + , + additive_combine_engine( + rnd1.engine(), + rnd2.engine(), + private_constructor_tag() + ) + ); + } + public: + SPROUT_CONSTEXPR additive_combine_engine() + : mlcg1_() + , mlcg2_() + {} + SPROUT_CONSTEXPR explicit additive_combine_engine(result_type const& seed) + : mlcg1_(seed) + , mlcg2_(seed) + {} + SPROUT_CONSTEXPR additive_combine_engine(typename MLCG1::result_type seed1, typename MLCG2::result_type seed2) + : mlcg1_(seed1) + , mlcg2_(seed2) + {} + SPROUT_CONSTEXPR result_type min() const { + return 1; + } + SPROUT_CONSTEXPR result_type max() const { + return first_base::modulus - 1; + } + SPROUT_CONSTEXPR sprout::random::random_result operator()() const { + return generate(mlcg1_(), mlcg2_()); + } + friend SPROUT_CONSTEXPR bool operator==(additive_combine_engine const& lhs, additive_combine_engine const& rhs) { + return lhs.mlcg1_ == rhs.mlcg1_ && lhs.mlcg2_ == rhs.mlcg2_; + } + friend SPROUT_CONSTEXPR bool operator!=(additive_combine_engine const& lhs, additive_combine_engine const& rhs) { + return !(lhs == rhs); + } + template + friend std::basic_istream& operator>>( + std::basic_istream& lhs, + additive_combine_engine& rhs + ) + { + return lhs >> rhs.mlcg1_ >> std::ws >> rhs.mlcg2_; + } + template + friend std::basic_ostream& operator<<( + std::basic_ostream& lhs, + additive_combine_engine const& rhs + ) + { + return lhs << rhs.mlcg1_ << ' ' << rhs.mlcg2_; + } + }; + + // + // ecuyer1988 + // + typedef additive_combine_engine< + sprout::random::linear_congruential_engine, + sprout::random::linear_congruential_engine + > ecuyer1988; + } // namespace random + + using sprout::random::ecuyer1988; +} // namespace sprout + +#endif // #ifndef SPROUT_RANDOM_ADDITIVE_COMBINE_HPP diff --git a/sprout/random/inversive_congruential.hpp b/sprout/random/inversive_congruential.hpp new file mode 100644 index 00000000..139231bb --- /dev/null +++ b/sprout/random/inversive_congruential.hpp @@ -0,0 +1,130 @@ +#ifndef SPROUT_RANDOM_INVERSIVE_CONGRUENTIAL_HPP +#define SPROUT_RANDOM_INVERSIVE_CONGRUENTIAL_HPP + +#include +#include +#include +#include +#include +#include + +namespace sprout { + namespace random { + // + // inversive_congruential_engine + // + template + class inversive_congruential_engine { + public: + typedef IntType result_type; + private: + struct private_constructor_tag {}; + public: + SPROUT_STATIC_CONSTEXPR IntType multiplier = a; + SPROUT_STATIC_CONSTEXPR IntType increment = b; + SPROUT_STATIC_CONSTEXPR IntType modulus = p; + SPROUT_STATIC_CONSTEXPR IntType default_seed = 1; + private: + static SPROUT_CONSTEXPR result_type static_min() { + return b == 0 ? 1 : 0; + } + static SPROUT_CONSTEXPR result_type static_max() { + return modulus - 1; + } + static SPROUT_CONSTEXPR bool arg_check_nothrow(IntType const& x0) { + return x0 >= static_min() && x0 <= static_max(); + } + static SPROUT_CONSTEXPR IntType arg_check(IntType const& x0) { + return arg_check_nothrow(x0) + ? x0 + : throw "assert(x0 >= static_min() && x0 <= static_max())" + ; + } + static SPROUT_CONSTEXPR IntType init_seed_2(IntType const& x0) { + return arg_check(increment == 0 && x0 == 0 ? 1 : x0); + } + static SPROUT_CONSTEXPR IntType init_seed_1(IntType const& x0) { + return init_seed_2(x0 <= 0 && x0 != 0 ? x0 + modulus : x0); + } + static SPROUT_CONSTEXPR IntType init_seed(IntType const& x0) { + return init_seed_1(modulus == 0 ? x0 : x0 % modulus); + } + private: + IntType x_; + private: + SPROUT_CONSTEXPR inversive_congruential_engine(IntType const& x, private_constructor_tag) + : x_(x) + {} + SPROUT_CONSTEXPR sprout::random::random_result generate(result_type result) const { + return sprout::random::random_result( + result, + inversive_congruential_engine(result, private_constructor_tag()) + ); + } + public: + SPROUT_CONSTEXPR inversive_congruential_engine() + : x_(init_seed(default_seed)) + {} + SPROUT_CONSTEXPR explicit inversive_congruential_engine(IntType const& x0) + : x_(init_seed(x0)) + {} + SPROUT_CONSTEXPR result_type min() const { + return static_min(); + } + SPROUT_CONSTEXPR result_type max() const { + return static_max(); + } + SPROUT_CONSTEXPR sprout::random::random_result operator()() const { + typedef sprout::random::detail::const_mod do_mod; + return generate(do_mod::mult_add(a, do_mod::invert(x_), b)); + } + friend SPROUT_CONSTEXPR bool operator==(inversive_congruential_engine const& lhs, inversive_congruential_engine const& rhs) { + return lhs.x_ == rhs.x_; + } + friend SPROUT_CONSTEXPR bool operator!=(inversive_congruential_engine const& lhs, inversive_congruential_engine const& rhs) { + return !(lhs == rhs); + } + template + friend std::basic_istream& operator>>( + std::basic_istream& lhs, + inversive_congruential_engine& rhs + ) + { + IntType x; + if (lhs >> x) { + if(arg_check_nothrow(x)) { + rhs.x_ = x; + } else { + lhs.setstate(std::ios_base::failbit); + } + } + return lhs; + } + template + friend std::basic_ostream& operator<<( + std::basic_ostream& lhs, + inversive_congruential_engine const& rhs + ) + { + return lhs << rhs.x_; + } + }; + template + SPROUT_CONSTEXPR IntType sprout::random::inversive_congruential_engine::multiplier; + template + SPROUT_CONSTEXPR IntType sprout::random::inversive_congruential_engine::increment; + template + SPROUT_CONSTEXPR IntType sprout::random::inversive_congruential_engine::modulus; + template + SPROUT_CONSTEXPR IntType sprout::random::inversive_congruential_engine::default_seed; + + // + // hellekalek1995 + // + typedef sprout::random::inversive_congruential_engine hellekalek1995; + } // namespace random + + using sprout::random::hellekalek1995; +} // namespace sprout + +#endif // #ifndef SPROUT_RANDOM_INVERSIVE_CONGRUENTIAL_HPP diff --git a/sprout/random/linear_congruential.hpp b/sprout/random/linear_congruential.hpp index 1b9cdb9a..96dc7bcc 100644 --- a/sprout/random/linear_congruential.hpp +++ b/sprout/random/linear_congruential.hpp @@ -29,6 +29,12 @@ namespace sprout { static_assert(m == 0 || a < m, "m == 0 || a < m"); static_assert(m == 0 || c < m, "m == 0 || c < m"); private: + static SPROUT_CONSTEXPR result_type static_min() { + return c == 0 ? 1 : 0; + } + static SPROUT_CONSTEXPR result_type static_max() { + return modulus - 1; + } static SPROUT_CONSTEXPR bool arg_check_nothrow(IntType const& x0) { return x0 >= static_min() && x0 <= static_max(); } @@ -47,13 +53,6 @@ namespace sprout { static SPROUT_CONSTEXPR IntType init_seed(IntType const& x0) { return init_seed_1(modulus == 0 ? x0 : x0 % modulus); } - public: - static SPROUT_CONSTEXPR result_type static_min() { - return c == 0 ? 1 : 0; - } - static SPROUT_CONSTEXPR result_type static_max() { - return modulus - 1; - } private: IntType x_; private: diff --git a/sprout/random/shuffle_order.hpp b/sprout/random/shuffle_order.hpp index 331a2fc5..55374e91 100644 --- a/sprout/random/shuffle_order.hpp +++ b/sprout/random/shuffle_order.hpp @@ -4,7 +4,6 @@ #include #include #include -#include #include #include #include @@ -141,20 +140,11 @@ namespace sprout { shuffle_order_engine& rhs ) { - using std::swap; - base_type rng; - sprout::array v; - result_type y; - lhs >> rng; + lhs >> rhs.rng_; for (std::size_t i = 0; i < k; ++i) { - lhs >> std::ws >> v[i]; - } - lhs >> std::ws >> y; - if (lhs) { - swap(rhs.rng_, rng); - swap(rhs.v_, v); - swap(rhs.y_, y); + lhs >> std::ws >> rhs.v_[i]; } + lhs >> std::ws >> rhs.y_; return lhs; } template