williamr@2: /* boost random/lagged_fibonacci.hpp header file williamr@2: * williamr@2: * Copyright Jens Maurer 2000-2001 williamr@2: * Distributed under the Boost Software License, Version 1.0. (See williamr@2: * accompanying file LICENSE_1_0.txt or copy at williamr@2: * http://www.boost.org/LICENSE_1_0.txt) williamr@2: * williamr@2: * See http://www.boost.org for most recent version including documentation. williamr@2: * williamr@2: * $Id: lagged_fibonacci.hpp,v 1.28 2005/05/21 15:57:00 dgregor Exp $ williamr@2: * williamr@2: * Revision history williamr@2: * 2001-02-18 moved to individual header files williamr@2: */ williamr@2: williamr@2: #ifndef BOOST_RANDOM_LAGGED_FIBONACCI_HPP williamr@2: #define BOOST_RANDOM_LAGGED_FIBONACCI_HPP williamr@2: williamr@2: #include williamr@2: #include williamr@2: #include // std::max williamr@2: #include williamr@2: #include // std::pow williamr@2: #include williamr@2: #include williamr@2: #include williamr@2: #include williamr@2: #include williamr@2: #include williamr@2: #include williamr@2: williamr@2: namespace boost { williamr@2: namespace random { williamr@2: williamr@2: #if BOOST_WORKAROUND(_MSC_FULL_VER, BOOST_TESTED_AT(13102292)) && BOOST_MSVC > 1300 williamr@2: # define BOOST_RANDOM_EXTRACT_LF williamr@2: #endif williamr@2: williamr@2: #if defined(__APPLE_CC__) && defined(__GNUC__) && (__GNUC__ == 3) && (__GNUC_MINOR__ <= 3) williamr@2: # define BOOST_RANDOM_EXTRACT_LF williamr@2: #endif williamr@2: williamr@2: # ifdef BOOST_RANDOM_EXTRACT_LF williamr@2: namespace detail williamr@2: { williamr@2: template williamr@2: IStream& williamr@2: extract_lagged_fibonacci_01( williamr@2: IStream& is williamr@2: , F const& f williamr@2: , unsigned int& i williamr@2: , RealType* x williamr@2: , RealType modulus) williamr@2: { williamr@2: is >> i >> std::ws; williamr@2: for(unsigned int i = 0; i < f.long_lag; ++i) williamr@2: { williamr@2: RealType value; williamr@2: is >> value >> std::ws; williamr@2: x[i] = value / modulus; williamr@2: } williamr@2: return is; williamr@2: } williamr@2: williamr@2: template williamr@2: IStream& williamr@2: extract_lagged_fibonacci( williamr@2: IStream& is williamr@2: , F const& f williamr@2: , unsigned int& i williamr@2: , UIntType* x) williamr@2: { williamr@2: is >> i >> std::ws; williamr@2: for(unsigned int i = 0; i < f.long_lag; ++i) williamr@2: is >> x[i] >> std::ws; williamr@2: return is; williamr@2: } williamr@2: } williamr@2: # endif williamr@2: williamr@2: template williamr@2: class lagged_fibonacci williamr@2: { williamr@2: public: williamr@2: typedef UIntType result_type; williamr@2: BOOST_STATIC_CONSTANT(bool, has_fixed_range = false); williamr@2: BOOST_STATIC_CONSTANT(int, word_size = w); williamr@2: BOOST_STATIC_CONSTANT(unsigned int, long_lag = p); williamr@2: BOOST_STATIC_CONSTANT(unsigned int, short_lag = q); williamr@2: williamr@2: result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; } williamr@2: result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return wordmask; } williamr@2: williamr@2: lagged_fibonacci() { init_wordmask(); seed(); } williamr@2: explicit lagged_fibonacci(uint32_t value) { init_wordmask(); seed(value); } williamr@2: template lagged_fibonacci(It& first, It last) williamr@2: { init_wordmask(); seed(first, last); } williamr@2: // compiler-generated copy ctor and assignment operator are fine williamr@2: williamr@2: private: williamr@2: void init_wordmask() williamr@2: { williamr@2: wordmask = 0; williamr@2: for(int i = 0; i < w; ++i) williamr@2: wordmask |= (1u << i); williamr@2: } williamr@2: williamr@2: public: williamr@2: void seed(uint32_t value = 331u) williamr@2: { williamr@2: minstd_rand0 gen(value); williamr@2: for(unsigned int j = 0; j < long_lag; ++j) williamr@2: x[j] = gen() & wordmask; williamr@2: i = long_lag; williamr@2: } williamr@2: williamr@2: template williamr@2: void seed(It& first, It last) williamr@2: { williamr@2: // word size could be smaller than the seed values williamr@2: unsigned int j; williamr@2: for(j = 0; j < long_lag && first != last; ++j, ++first) williamr@2: x[j] = *first & wordmask; williamr@2: i = long_lag; williamr@2: if(first == last && j < long_lag) williamr@2: throw std::invalid_argument("lagged_fibonacci::seed"); williamr@2: } williamr@2: williamr@2: result_type operator()() williamr@2: { williamr@2: if(i >= long_lag) williamr@2: fill(); williamr@2: return x[i++]; williamr@2: } williamr@2: williamr@2: static bool validation(result_type x) williamr@2: { williamr@2: return x == val; williamr@2: } williamr@2: williamr@2: #ifndef BOOST_NO_OPERATORS_IN_NAMESPACE williamr@2: williamr@2: #ifndef BOOST_NO_MEMBER_TEMPLATE_FRIENDS williamr@2: template williamr@2: friend std::basic_ostream& williamr@2: operator<<(std::basic_ostream& os, const lagged_fibonacci& f) williamr@2: { williamr@2: os << f.i << " "; williamr@2: for(unsigned int i = 0; i < f.long_lag; ++i) williamr@2: os << f.x[i] << " "; williamr@2: return os; williamr@2: } williamr@2: williamr@2: template williamr@2: friend std::basic_istream& williamr@2: operator>>(std::basic_istream& is, lagged_fibonacci& f) williamr@2: { williamr@2: # ifdef BOOST_RANDOM_EXTRACT_LF williamr@2: return detail::extract_lagged_fibonacci(is, f, f.i, f.x); williamr@2: # else williamr@2: is >> f.i >> std::ws; williamr@2: for(unsigned int i = 0; i < f.long_lag; ++i) williamr@2: is >> f.x[i] >> std::ws; williamr@2: return is; williamr@2: # endif williamr@2: } williamr@2: #endif williamr@2: williamr@2: friend bool operator==(const lagged_fibonacci& x, const lagged_fibonacci& y) williamr@2: { return x.i == y.i && std::equal(x.x, x.x+long_lag, y.x); } williamr@2: friend bool operator!=(const lagged_fibonacci& x, williamr@2: const lagged_fibonacci& y) williamr@2: { return !(x == y); } williamr@2: #else williamr@2: // Use a member function; Streamable concept not supported. williamr@2: bool operator==(const lagged_fibonacci& rhs) const williamr@2: { return i == rhs.i && std::equal(x, x+long_lag, rhs.x); } williamr@2: bool operator!=(const lagged_fibonacci& rhs) const williamr@2: { return !(*this == rhs); } williamr@2: #endif williamr@2: williamr@2: private: williamr@2: void fill(); williamr@2: UIntType wordmask; williamr@2: unsigned int i; williamr@2: UIntType x[long_lag]; williamr@2: }; williamr@2: williamr@2: #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION williamr@2: // A definition is required even for integral static constants williamr@2: template williamr@2: const bool lagged_fibonacci::has_fixed_range; williamr@2: template williamr@2: const unsigned int lagged_fibonacci::long_lag; williamr@2: template williamr@2: const unsigned int lagged_fibonacci::short_lag; williamr@2: #endif williamr@2: williamr@2: template williamr@2: void lagged_fibonacci::fill() williamr@2: { williamr@2: // two loops to avoid costly modulo operations williamr@2: { // extra scope for MSVC brokenness w.r.t. for scope williamr@2: for(unsigned int j = 0; j < short_lag; ++j) williamr@2: x[j] = (x[j] + x[j+(long_lag-short_lag)]) & wordmask; williamr@2: } williamr@2: for(unsigned int j = short_lag; j < long_lag; ++j) williamr@2: x[j] = (x[j] + x[j-short_lag]) & wordmask; williamr@2: i = 0; williamr@2: } williamr@2: williamr@2: williamr@2: williamr@2: // lagged Fibonacci generator for the range [0..1) williamr@2: // contributed by Matthias Troyer williamr@2: // for p=55, q=24 originally by G. J. Mitchell and D. P. Moore 1958 williamr@2: williamr@2: template williamr@2: struct fibonacci_validation williamr@2: { williamr@2: BOOST_STATIC_CONSTANT(bool, is_specialized = false); williamr@2: static T value() { return 0; } williamr@2: static T tolerance() { return 0; } williamr@2: }; williamr@2: williamr@2: #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION williamr@2: // A definition is required even for integral static constants williamr@2: template williamr@2: const bool fibonacci_validation::is_specialized; williamr@2: #endif williamr@2: williamr@2: #define BOOST_RANDOM_FIBONACCI_VAL(T,P,Q,V,E) \ williamr@2: template<> \ williamr@2: struct fibonacci_validation \ williamr@2: { \ williamr@2: BOOST_STATIC_CONSTANT(bool, is_specialized = true); \ williamr@2: static T value() { return V; } \ williamr@2: static T tolerance() \ williamr@2: { return (std::max)(E, static_cast(5*std::numeric_limits::epsilon())); } \ williamr@2: }; williamr@2: // (The extra static_cast in the std::max call above is actually williamr@2: // unnecessary except for HP aCC 1.30, which claims that williamr@2: // numeric_limits::epsilon() doesn't actually return a double.) williamr@2: williamr@2: BOOST_RANDOM_FIBONACCI_VAL(double, 607, 273, 0.4293817707235914, 1e-14) williamr@2: BOOST_RANDOM_FIBONACCI_VAL(double, 1279, 418, 0.9421630240437659, 1e-14) williamr@2: BOOST_RANDOM_FIBONACCI_VAL(double, 2281, 1252, 0.1768114046909004, 1e-14) williamr@2: BOOST_RANDOM_FIBONACCI_VAL(double, 3217, 576, 0.1956232694868209, 1e-14) williamr@2: BOOST_RANDOM_FIBONACCI_VAL(double, 4423, 2098, 0.9499762202147172, 1e-14) williamr@2: BOOST_RANDOM_FIBONACCI_VAL(double, 9689, 5502, 0.05737836943695162, 1e-14) williamr@2: BOOST_RANDOM_FIBONACCI_VAL(double, 19937, 9842, 0.5076528587449834, 1e-14) williamr@2: BOOST_RANDOM_FIBONACCI_VAL(double, 23209, 13470, 0.5414473810619185, 1e-14) williamr@2: BOOST_RANDOM_FIBONACCI_VAL(double, 44497,21034, 0.254135073399297, 1e-14) williamr@2: williamr@2: #undef BOOST_RANDOM_FIBONACCI_VAL williamr@2: williamr@2: template williamr@2: class lagged_fibonacci_01 williamr@2: { williamr@2: public: williamr@2: typedef RealType result_type; williamr@2: BOOST_STATIC_CONSTANT(bool, has_fixed_range = false); williamr@2: BOOST_STATIC_CONSTANT(int, word_size = w); williamr@2: BOOST_STATIC_CONSTANT(unsigned int, long_lag = p); williamr@2: BOOST_STATIC_CONSTANT(unsigned int, short_lag = q); williamr@2: williamr@2: lagged_fibonacci_01() { init_modulus(); seed(); } williamr@2: explicit lagged_fibonacci_01(uint32_t value) { init_modulus(); seed(value); } williamr@2: template williamr@2: explicit lagged_fibonacci_01(Generator & gen) { init_modulus(); seed(gen); } williamr@2: template lagged_fibonacci_01(It& first, It last) williamr@2: { init_modulus(); seed(first, last); } williamr@2: // compiler-generated copy ctor and assignment operator are fine williamr@2: williamr@2: private: williamr@2: void init_modulus() williamr@2: { williamr@2: #ifndef BOOST_NO_STDC_NAMESPACE williamr@2: // allow for Koenig lookup williamr@2: using std::pow; williamr@2: #endif williamr@2: _modulus = pow(RealType(2), word_size); williamr@2: } williamr@2: williamr@2: public: williamr@2: void seed(uint32_t value = 331u) williamr@2: { williamr@2: minstd_rand0 intgen(value); williamr@2: seed(intgen); williamr@2: } williamr@2: williamr@2: // For GCC, moving this function out-of-line prevents inlining, which may williamr@2: // reduce overall object code size. However, MSVC does not grok williamr@2: // out-of-line template member functions. williamr@2: template williamr@2: void seed(Generator & gen) williamr@2: { williamr@2: // use pass-by-reference, but wrap argument in pass_through_engine williamr@2: typedef detail::pass_through_engine ref_gen; williamr@2: uniform_01 gen01 = williamr@2: uniform_01(ref_gen(gen)); williamr@2: // I could have used std::generate_n, but it takes "gen" by value williamr@2: for(unsigned int j = 0; j < long_lag; ++j) williamr@2: x[j] = gen01(); williamr@2: i = long_lag; williamr@2: } williamr@2: williamr@2: template williamr@2: void seed(It& first, It last) williamr@2: { williamr@2: #ifndef BOOST_NO_STDC_NAMESPACE williamr@2: // allow for Koenig lookup williamr@2: using std::fmod; williamr@2: using std::pow; williamr@2: #endif williamr@2: unsigned long mask = ~((~0u) << (w%32)); // now lowest w bits set williamr@2: RealType two32 = pow(RealType(2), 32); williamr@2: unsigned int j; williamr@2: for(j = 0; j < long_lag && first != last; ++j, ++first) { williamr@2: x[j] = RealType(0); williamr@2: for(int i = 0; i < w/32 && first != last; ++i, ++first) williamr@2: x[j] += *first / pow(two32,i+1); williamr@2: if(first != last && mask != 0) williamr@2: x[j] += fmod((*first & mask) / _modulus, RealType(1)); williamr@2: } williamr@2: i = long_lag; williamr@2: if(first == last && j < long_lag) williamr@2: throw std::invalid_argument("lagged_fibonacci_01::seed"); williamr@2: } williamr@2: williamr@2: result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(0); } williamr@2: result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(1); } williamr@2: williamr@2: result_type operator()() williamr@2: { williamr@2: if(i >= long_lag) williamr@2: fill(); williamr@2: return x[i++]; williamr@2: } williamr@2: williamr@2: static bool validation(result_type x) williamr@2: { williamr@2: result_type v = fibonacci_validation::value(); williamr@2: result_type epsilon = fibonacci_validation::tolerance(); williamr@2: // std::abs is a source of trouble: sometimes, it's not overloaded williamr@2: // for double, plus the usual namespace std noncompliance -> avoid it williamr@2: // using std::abs; williamr@2: // return abs(x - v) < 5 * epsilon williamr@2: return x > v - epsilon && x < v + epsilon; williamr@2: } williamr@2: williamr@2: #ifndef BOOST_NO_OPERATORS_IN_NAMESPACE williamr@2: williamr@2: #ifndef BOOST_NO_MEMBER_TEMPLATE_FRIENDS williamr@2: template williamr@2: friend std::basic_ostream& williamr@2: operator<<(std::basic_ostream& os, const lagged_fibonacci_01&f) williamr@2: { williamr@2: #ifndef BOOST_NO_STDC_NAMESPACE williamr@2: // allow for Koenig lookup williamr@2: using std::pow; williamr@2: #endif williamr@2: os << f.i << " "; williamr@2: std::ios_base::fmtflags oldflags = os.flags(os.dec | os.fixed | os.left); williamr@2: for(unsigned int i = 0; i < f.long_lag; ++i) williamr@2: os << f.x[i] * f._modulus << " "; williamr@2: os.flags(oldflags); williamr@2: return os; williamr@2: } williamr@2: williamr@2: template williamr@2: friend std::basic_istream& williamr@2: operator>>(std::basic_istream& is, lagged_fibonacci_01& f) williamr@2: { williamr@2: # ifdef BOOST_RANDOM_EXTRACT_LF williamr@2: return detail::extract_lagged_fibonacci_01(is, f, f.i, f.x, f._modulus); williamr@2: # else williamr@2: is >> f.i >> std::ws; williamr@2: for(unsigned int i = 0; i < f.long_lag; ++i) { williamr@2: typename lagged_fibonacci_01::result_type value; williamr@2: is >> value >> std::ws; williamr@2: f.x[i] = value / f._modulus; williamr@2: } williamr@2: return is; williamr@2: # endif williamr@2: } williamr@2: #endif williamr@2: williamr@2: friend bool operator==(const lagged_fibonacci_01& x, williamr@2: const lagged_fibonacci_01& y) williamr@2: { return x.i == y.i && std::equal(x.x, x.x+long_lag, y.x); } williamr@2: friend bool operator!=(const lagged_fibonacci_01& x, williamr@2: const lagged_fibonacci_01& y) williamr@2: { return !(x == y); } williamr@2: #else williamr@2: // Use a member function; Streamable concept not supported. williamr@2: bool operator==(const lagged_fibonacci_01& rhs) const williamr@2: { return i == rhs.i && std::equal(x, x+long_lag, rhs.x); } williamr@2: bool operator!=(const lagged_fibonacci_01& rhs) const williamr@2: { return !(*this == rhs); } williamr@2: #endif williamr@2: williamr@2: private: williamr@2: void fill(); williamr@2: unsigned int i; williamr@2: RealType x[long_lag]; williamr@2: RealType _modulus; williamr@2: }; williamr@2: williamr@2: #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION williamr@2: // A definition is required even for integral static constants williamr@2: template williamr@2: const bool lagged_fibonacci_01::has_fixed_range; williamr@2: template williamr@2: const unsigned int lagged_fibonacci_01::long_lag; williamr@2: template williamr@2: const unsigned int lagged_fibonacci_01::short_lag; williamr@2: template williamr@2: const int lagged_fibonacci_01::word_size; williamr@2: williamr@2: #endif williamr@2: williamr@2: template williamr@2: void lagged_fibonacci_01::fill() williamr@2: { williamr@2: // two loops to avoid costly modulo operations williamr@2: { // extra scope for MSVC brokenness w.r.t. for scope williamr@2: for(unsigned int j = 0; j < short_lag; ++j) { williamr@2: RealType t = x[j] + x[j+(long_lag-short_lag)]; williamr@2: if(t >= RealType(1)) williamr@2: t -= RealType(1); williamr@2: x[j] = t; williamr@2: } williamr@2: } williamr@2: for(unsigned int j = short_lag; j < long_lag; ++j) { williamr@2: RealType t = x[j] + x[j-short_lag]; williamr@2: if(t >= RealType(1)) williamr@2: t -= RealType(1); williamr@2: x[j] = t; williamr@2: } williamr@2: i = 0; williamr@2: } williamr@2: williamr@2: } // namespace random williamr@2: williamr@2: typedef random::lagged_fibonacci_01 lagged_fibonacci607; williamr@2: typedef random::lagged_fibonacci_01 lagged_fibonacci1279; williamr@2: typedef random::lagged_fibonacci_01 lagged_fibonacci2281; williamr@2: typedef random::lagged_fibonacci_01 lagged_fibonacci3217; williamr@2: typedef random::lagged_fibonacci_01 lagged_fibonacci4423; williamr@2: typedef random::lagged_fibonacci_01 lagged_fibonacci9689; williamr@2: typedef random::lagged_fibonacci_01 lagged_fibonacci19937; williamr@2: typedef random::lagged_fibonacci_01 lagged_fibonacci23209; williamr@2: typedef random::lagged_fibonacci_01 lagged_fibonacci44497; williamr@2: williamr@2: williamr@2: // It is possible to partially specialize uniform_01<> on lagged_fibonacci_01<> williamr@2: // to help the compiler generate efficient code. For GCC, this seems useless, williamr@2: // because GCC optimizes (x-0)/(1-0) to (x-0). This is good enough for now. williamr@2: williamr@2: } // namespace boost williamr@2: williamr@2: #endif // BOOST_RANDOM_LAGGED_FIBONACCI_HPP