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