| /* boost random/subtract_with_carry.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: subtract_with_carry.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ |
| * |
| * Revision history |
| * 2002-03-02 created |
| */ |
| |
| #ifndef BOOST_RANDOM_SUBTRACT_WITH_CARRY_HPP |
| #define BOOST_RANDOM_SUBTRACT_WITH_CARRY_HPP |
| |
| #include <boost/config/no_tr1/cmath.hpp> |
| #include <iostream> |
| #include <algorithm> // std::equal |
| #include <stdexcept> |
| #include <boost/config/no_tr1/cmath.hpp> // std::pow |
| #include <boost/config.hpp> |
| #include <boost/limits.hpp> |
| #include <boost/cstdint.hpp> |
| #include <boost/static_assert.hpp> |
| #include <boost/detail/workaround.hpp> |
| #include <boost/random/detail/config.hpp> |
| #include <boost/random/detail/seed.hpp> |
| #include <boost/random/linear_congruential.hpp> |
| |
| |
| namespace boost { |
| namespace random { |
| |
| #if BOOST_WORKAROUND(_MSC_FULL_VER, BOOST_TESTED_AT(13102292)) && BOOST_MSVC > 1300 |
| # define BOOST_RANDOM_EXTRACT_SWC_01 |
| #endif |
| |
| #if defined(__APPLE_CC__) && defined(__GNUC__) && (__GNUC__ == 3) && (__GNUC_MINOR__ <= 3) |
| # define BOOST_RANDOM_EXTRACT_SWC_01 |
| #endif |
| |
| # ifdef BOOST_RANDOM_EXTRACT_SWC_01 |
| namespace detail |
| { |
| template <class IStream, class SubtractWithCarry, class RealType> |
| void extract_subtract_with_carry_01( |
| IStream& is |
| , SubtractWithCarry& f |
| , RealType& carry |
| , RealType* x |
| , RealType modulus) |
| { |
| RealType value; |
| for(unsigned int j = 0; j < f.long_lag; ++j) { |
| is >> value >> std::ws; |
| x[j] = value / modulus; |
| } |
| is >> value >> std::ws; |
| carry = value / modulus; |
| } |
| } |
| # endif |
| |
| /** |
| * Instantiations of @c subtract_with_carry model a |
| * \pseudo_random_number_generator. The algorithm is |
| * described in |
| * |
| * @blockquote |
| * "A New Class of Random Number Generators", George |
| * Marsaglia and Arif Zaman, Annals of Applied Probability, |
| * Volume 1, Number 3 (1991), 462-480. |
| * @endblockquote |
| */ |
| template<class IntType, IntType m, unsigned int s, unsigned int r, |
| IntType val> |
| class subtract_with_carry |
| { |
| public: |
| typedef IntType result_type; |
| BOOST_STATIC_CONSTANT(bool, has_fixed_range = true); |
| BOOST_STATIC_CONSTANT(result_type, min_value = 0); |
| BOOST_STATIC_CONSTANT(result_type, max_value = m-1); |
| BOOST_STATIC_CONSTANT(result_type, modulus = m); |
| BOOST_STATIC_CONSTANT(unsigned int, long_lag = r); |
| BOOST_STATIC_CONSTANT(unsigned int, short_lag = s); |
| |
| subtract_with_carry() { |
| // MSVC fails BOOST_STATIC_ASSERT with std::numeric_limits at class scope |
| #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS |
| BOOST_STATIC_ASSERT(std::numeric_limits<result_type>::is_signed); |
| BOOST_STATIC_ASSERT(std::numeric_limits<result_type>::is_integer); |
| #endif |
| seed(); |
| } |
| BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(subtract_with_carry, uint32_t, value) |
| { seed(value); } |
| BOOST_RANDOM_DETAIL_GENERATOR_CONSTRUCTOR(subtract_with_carry, Generator, gen) |
| { seed(gen); } |
| template<class It> subtract_with_carry(It& first, It last) { seed(first,last); } |
| |
| // compiler-generated copy ctor and assignment operator are fine |
| |
| void seed() { seed(19780503u); } |
| BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(subtract_with_carry, uint32_t, value) |
| { |
| random::linear_congruential<int32_t, 40014, 0, 2147483563, 0> intgen(value); |
| seed(intgen); |
| } |
| |
| // For GCC, moving this function out-of-line prevents inlining, which may |
| // reduce overall object code size. However, MSVC does not grok |
| // out-of-line template member functions. |
| BOOST_RANDOM_DETAIL_GENERATOR_SEED(subtract_with_carry, Generator, gen) |
| { |
| // I could have used std::generate_n, but it takes "gen" by value |
| for(unsigned int j = 0; j < long_lag; ++j) |
| x[j] = gen() % modulus; |
| carry = (x[long_lag-1] == 0); |
| k = 0; |
| } |
| |
| template<class It> |
| void seed(It& first, It last) |
| { |
| unsigned int j; |
| for(j = 0; j < long_lag && first != last; ++j, ++first) |
| x[j] = *first % modulus; |
| if(first == last && j < long_lag) |
| throw std::invalid_argument("subtract_with_carry::seed"); |
| carry = (x[long_lag-1] == 0); |
| k = 0; |
| } |
| |
| result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return min_value; } |
| result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return max_value; } |
| |
| result_type operator()() |
| { |
| int short_index = k - short_lag; |
| if(short_index < 0) |
| short_index += long_lag; |
| IntType delta; |
| if (x[short_index] >= x[k] + carry) { |
| // x(n) >= 0 |
| delta = x[short_index] - (x[k] + carry); |
| carry = 0; |
| } else { |
| // x(n) < 0 |
| delta = modulus - x[k] - carry + x[short_index]; |
| carry = 1; |
| } |
| x[k] = delta; |
| ++k; |
| if(k >= long_lag) |
| k = 0; |
| return delta; |
| } |
| |
| public: |
| static bool validation(result_type x) { return x == val; } |
| |
| #ifndef BOOST_NO_OPERATORS_IN_NAMESPACE |
| |
| #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 subtract_with_carry& f) |
| { |
| for(unsigned int j = 0; j < f.long_lag; ++j) |
| os << f.compute(j) << " "; |
| os << f.carry << " "; |
| return os; |
| } |
| |
| template<class CharT, class Traits> |
| friend std::basic_istream<CharT,Traits>& |
| operator>>(std::basic_istream<CharT,Traits>& is, subtract_with_carry& f) |
| { |
| for(unsigned int j = 0; j < f.long_lag; ++j) |
| is >> f.x[j] >> std::ws; |
| is >> f.carry >> std::ws; |
| f.k = 0; |
| return is; |
| } |
| #endif |
| |
| friend bool operator==(const subtract_with_carry& x, const subtract_with_carry& y) |
| { |
| for(unsigned int j = 0; j < r; ++j) |
| if(x.compute(j) != y.compute(j)) |
| return false; |
| return true; |
| } |
| |
| friend bool operator!=(const subtract_with_carry& x, const subtract_with_carry& y) |
| { return !(x == y); } |
| #else |
| // Use a member function; Streamable concept not supported. |
| bool operator==(const subtract_with_carry& rhs) const |
| { |
| for(unsigned int j = 0; j < r; ++j) |
| if(compute(j) != rhs.compute(j)) |
| return false; |
| return true; |
| } |
| |
| bool operator!=(const subtract_with_carry& rhs) const |
| { return !(*this == rhs); } |
| #endif |
| |
| private: |
| /// \cond hide_private_members |
| // returns x(i-r+index), where index is in 0..r-1 |
| IntType compute(unsigned int index) const |
| { |
| return x[(k+index) % long_lag]; |
| } |
| /// \endcond |
| |
| // state representation; next output (state) is x(i) |
| // x[0] ... x[k] x[k+1] ... x[long_lag-1] represents |
| // x(i-k) ... x(i) x(i+1) ... x(i-k+long_lag-1) |
| // speed: base: 20-25 nsec |
| // ranlux_4: 230 nsec, ranlux_7: 430 nsec, ranlux_14: 810 nsec |
| // This state representation makes operator== and save/restore more |
| // difficult, because we've already computed "too much" and thus |
| // have to undo some steps to get at x(i-r) etc. |
| |
| // state representation: next output (state) is x(i) |
| // x[0] ... x[k] x[k+1] ... x[long_lag-1] represents |
| // x(i-k) ... x(i) x(i-long_lag+1) ... x(i-k-1) |
| // speed: base 28 nsec |
| // ranlux_4: 370 nsec, ranlux_7: 688 nsec, ranlux_14: 1343 nsec |
| IntType x[long_lag]; |
| unsigned int k; |
| int carry; |
| }; |
| |
| #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION |
| // A definition is required even for integral static constants |
| template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val> |
| const bool subtract_with_carry<IntType, m, s, r, val>::has_fixed_range; |
| template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val> |
| const IntType subtract_with_carry<IntType, m, s, r, val>::min_value; |
| template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val> |
| const IntType subtract_with_carry<IntType, m, s, r, val>::max_value; |
| template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val> |
| const IntType subtract_with_carry<IntType, m, s, r, val>::modulus; |
| template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val> |
| const unsigned int subtract_with_carry<IntType, m, s, r, val>::long_lag; |
| template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val> |
| const unsigned int subtract_with_carry<IntType, m, s, r, val>::short_lag; |
| #endif |
| |
| |
| // use a floating-point representation to produce values in [0..1) |
| /** @copydoc boost::random::subtract_with_carry */ |
| template<class RealType, int w, unsigned int s, unsigned int r, int val=0> |
| class subtract_with_carry_01 |
| { |
| public: |
| typedef RealType result_type; |
| BOOST_STATIC_CONSTANT(bool, has_fixed_range = false); |
| BOOST_STATIC_CONSTANT(int, word_size = w); |
| BOOST_STATIC_CONSTANT(unsigned int, long_lag = r); |
| BOOST_STATIC_CONSTANT(unsigned int, short_lag = s); |
| |
| #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS |
| BOOST_STATIC_ASSERT(!std::numeric_limits<result_type>::is_integer); |
| #endif |
| |
| subtract_with_carry_01() { init_modulus(); seed(); } |
| explicit subtract_with_carry_01(uint32_t value) |
| { init_modulus(); seed(value); } |
| template<class It> subtract_with_carry_01(It& first, It last) |
| { init_modulus(); seed(first,last); } |
| |
| private: |
| /// \cond hide_private_members |
| void init_modulus() |
| { |
| #ifndef BOOST_NO_STDC_NAMESPACE |
| // allow for Koenig lookup |
| using std::pow; |
| #endif |
| _modulus = pow(RealType(2), word_size); |
| } |
| /// \endcond hide_private_members |
| |
| public: |
| // compiler-generated copy ctor and assignment operator are fine |
| |
| void seed(uint32_t value = 19780503u) |
| { |
| #ifndef BOOST_NO_STDC_NAMESPACE |
| // allow for Koenig lookup |
| using std::fmod; |
| #endif |
| random::linear_congruential<int32_t, 40014, 0, 2147483563, 0> gen(value); |
| unsigned long array[(w+31)/32 * long_lag]; |
| for(unsigned int j = 0; j < sizeof(array)/sizeof(unsigned long); ++j) |
| array[j] = gen(); |
| unsigned long * start = array; |
| seed(start, array + sizeof(array)/sizeof(unsigned long)); |
| } |
| |
| template<class It> |
| void seed(It& first, It last) |
| { |
| #ifndef BOOST_NO_STDC_NAMESPACE |
| // allow for Koenig lookup |
| using std::fmod; |
| using std::pow; |
| #endif |
| unsigned long mask = ~((~0u) << (w%32)); // now lowest (w%32) bits set |
| RealType two32 = pow(RealType(2), 32); |
| unsigned int j; |
| for(j = 0; j < long_lag && first != last; ++j) { |
| x[j] = RealType(0); |
| for(int i = 0; i < w/32 && first != last; ++i, ++first) |
| x[j] += *first / pow(two32,i+1); |
| if(first != last && mask != 0) { |
| x[j] += fmod((*first & mask) / _modulus, RealType(1)); |
| ++first; |
| } |
| } |
| if(first == last && j < long_lag) |
| throw std::invalid_argument("subtract_with_carry_01::seed"); |
| carry = (x[long_lag-1] ? 0 : 1 / _modulus); |
| k = 0; |
| } |
| |
| result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(0); } |
| result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(1); } |
| |
| result_type operator()() |
| { |
| int short_index = k - short_lag; |
| if(short_index < 0) |
| short_index += long_lag; |
| RealType delta = x[short_index] - x[k] - carry; |
| if(delta < 0) { |
| delta += RealType(1); |
| carry = RealType(1)/_modulus; |
| } else { |
| carry = 0; |
| } |
| x[k] = delta; |
| ++k; |
| if(k >= long_lag) |
| k = 0; |
| return delta; |
| } |
| |
| static bool validation(result_type x) |
| { return x == val/pow(RealType(2), word_size); } |
| |
| #ifndef BOOST_NO_OPERATORS_IN_NAMESPACE |
| |
| #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 subtract_with_carry_01& f) |
| { |
| #ifndef BOOST_NO_STDC_NAMESPACE |
| // allow for Koenig lookup |
| using std::pow; |
| #endif |
| std::ios_base::fmtflags oldflags = os.flags(os.dec | os.fixed | os.left); |
| for(unsigned int j = 0; j < f.long_lag; ++j) |
| os << (f.compute(j) * f._modulus) << " "; |
| os << (f.carry * f._modulus); |
| os.flags(oldflags); |
| return os; |
| } |
| |
| template<class CharT, class Traits> |
| friend std::basic_istream<CharT,Traits>& |
| operator>>(std::basic_istream<CharT,Traits>& is, subtract_with_carry_01& f) |
| { |
| # ifdef BOOST_RANDOM_EXTRACT_SWC_01 |
| detail::extract_subtract_with_carry_01(is, f, f.carry, f.x, f._modulus); |
| # else |
| // MSVC (up to 7.1) and Borland (up to 5.64) don't handle the template type |
| // parameter "RealType" available from the class template scope, so use |
| // the member typedef |
| typename subtract_with_carry_01::result_type value; |
| for(unsigned int j = 0; j < long_lag; ++j) { |
| is >> value >> std::ws; |
| f.x[j] = value / f._modulus; |
| } |
| is >> value >> std::ws; |
| f.carry = value / f._modulus; |
| # endif |
| f.k = 0; |
| return is; |
| } |
| #endif |
| |
| friend bool operator==(const subtract_with_carry_01& x, |
| const subtract_with_carry_01& y) |
| { |
| for(unsigned int j = 0; j < r; ++j) |
| if(x.compute(j) != y.compute(j)) |
| return false; |
| return true; |
| } |
| |
| friend bool operator!=(const subtract_with_carry_01& x, |
| const subtract_with_carry_01& y) |
| { return !(x == y); } |
| #else |
| // Use a member function; Streamable concept not supported. |
| bool operator==(const subtract_with_carry_01& rhs) const |
| { |
| for(unsigned int j = 0; j < r; ++j) |
| if(compute(j) != rhs.compute(j)) |
| return false; |
| return true; |
| } |
| |
| bool operator!=(const subtract_with_carry_01& rhs) const |
| { return !(*this == rhs); } |
| #endif |
| |
| private: |
| /// \cond hide_private_members |
| RealType compute(unsigned int index) const; |
| /// \endcond |
| unsigned int k; |
| RealType carry; |
| RealType x[long_lag]; |
| RealType _modulus; |
| }; |
| |
| #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION |
| // A definition is required even for integral static constants |
| template<class RealType, int w, unsigned int s, unsigned int r, int val> |
| const bool subtract_with_carry_01<RealType, w, s, r, val>::has_fixed_range; |
| template<class RealType, int w, unsigned int s, unsigned int r, int val> |
| const int subtract_with_carry_01<RealType, w, s, r, val>::word_size; |
| template<class RealType, int w, unsigned int s, unsigned int r, int val> |
| const unsigned int subtract_with_carry_01<RealType, w, s, r, val>::long_lag; |
| template<class RealType, int w, unsigned int s, unsigned int r, int val> |
| const unsigned int subtract_with_carry_01<RealType, w, s, r, val>::short_lag; |
| #endif |
| |
| /// \cond hide_private_members |
| template<class RealType, int w, unsigned int s, unsigned int r, int val> |
| RealType subtract_with_carry_01<RealType, w, s, r, val>::compute(unsigned int index) const |
| { |
| return x[(k+index) % long_lag]; |
| } |
| /// \endcond |
| |
| } // namespace random |
| } // namespace boost |
| |
| #endif // BOOST_RANDOM_SUBTRACT_WITH_CARRY_HPP |