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 <cmath>
williamr@2: #include <iostream>
williamr@2: #include <algorithm>     // std::max
williamr@2: #include <iterator>
williamr@2: #include <cmath>         // std::pow
williamr@2: #include <boost/config.hpp>
williamr@2: #include <boost/limits.hpp>
williamr@2: #include <boost/cstdint.hpp>
williamr@2: #include <boost/detail/workaround.hpp>
williamr@2: #include <boost/random/linear_congruential.hpp>
williamr@2: #include <boost/random/uniform_01.hpp>
williamr@2: #include <boost/random/detail/pass_through_engine.hpp>
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<class IStream, class F, class RealType>
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<class IStream, class F, class UIntType>
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<class UIntType, int w, unsigned int p, unsigned int q,
williamr@2:          UIntType val = 0>
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<class It> 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<class It>
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<class CharT, class Traits>
williamr@2:   friend std::basic_ostream<CharT,Traits>&
williamr@2:   operator<<(std::basic_ostream<CharT,Traits>& 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<class CharT, class Traits>
williamr@2:   friend std::basic_istream<CharT, Traits>&
williamr@2:   operator>>(std::basic_istream<CharT, Traits>& 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<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
williamr@2: const bool lagged_fibonacci<UIntType, w, p, q, val>::has_fixed_range;
williamr@2: template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
williamr@2: const unsigned int lagged_fibonacci<UIntType, w, p, q, val>::long_lag;
williamr@2: template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
williamr@2: const unsigned int lagged_fibonacci<UIntType, w, p, q, val>::short_lag;
williamr@2: #endif
williamr@2: 
williamr@2: template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
williamr@2: void lagged_fibonacci<UIntType, w, p, q, val>::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<class T, unsigned int p, unsigned int q>
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<class T, unsigned int p, unsigned int q>
williamr@2: const bool fibonacci_validation<T, p, q>::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<T, P, Q>  \
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<T>(5*std::numeric_limits<T>::epsilon())); } \
williamr@2: };
williamr@2: // (The extra static_cast<T> in the std::max call above is actually
williamr@2: // unnecessary except for HP aCC 1.30, which claims that
williamr@2: // numeric_limits<double>::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<class RealType, int w, unsigned int p, unsigned int q>
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<class Generator>
williamr@2:   explicit lagged_fibonacci_01(Generator & gen) { init_modulus(); seed(gen); }
williamr@2:   template<class It> 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<class Generator>
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<Generator&> ref_gen;
williamr@2:     uniform_01<ref_gen, RealType> gen01 =
williamr@2:       uniform_01<ref_gen, RealType>(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<class It>
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<result_type, p, q>::value();
williamr@2:     result_type epsilon = fibonacci_validation<result_type, p, q>::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<class CharT, class Traits>
williamr@2:   friend std::basic_ostream<CharT,Traits>&
williamr@2:   operator<<(std::basic_ostream<CharT,Traits>& 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<class CharT, class Traits>
williamr@2:   friend std::basic_istream<CharT, Traits>&
williamr@2:   operator>>(std::basic_istream<CharT, Traits>& 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<class RealType, int w, unsigned int p, unsigned int q>
williamr@2: const bool lagged_fibonacci_01<RealType, w, p, q>::has_fixed_range;
williamr@2: template<class RealType, int w, unsigned int p, unsigned int q>
williamr@2: const unsigned int lagged_fibonacci_01<RealType, w, p, q>::long_lag;
williamr@2: template<class RealType, int w, unsigned int p, unsigned int q>
williamr@2: const unsigned int lagged_fibonacci_01<RealType, w, p, q>::short_lag;
williamr@2: template<class RealType, int w, unsigned int p, unsigned int q>
williamr@2: const int lagged_fibonacci_01<RealType,w,p,q>::word_size;
williamr@2: 
williamr@2: #endif
williamr@2: 
williamr@2: template<class RealType, int w, unsigned int p, unsigned int q>
williamr@2: void lagged_fibonacci_01<RealType, w, p, q>::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<double, 48, 607, 273> lagged_fibonacci607;
williamr@2: typedef random::lagged_fibonacci_01<double, 48, 1279, 418> lagged_fibonacci1279;
williamr@2: typedef random::lagged_fibonacci_01<double, 48, 2281, 1252> lagged_fibonacci2281;
williamr@2: typedef random::lagged_fibonacci_01<double, 48, 3217, 576> lagged_fibonacci3217;
williamr@2: typedef random::lagged_fibonacci_01<double, 48, 4423, 2098> lagged_fibonacci4423;
williamr@2: typedef random::lagged_fibonacci_01<double, 48, 9689, 5502> lagged_fibonacci9689;
williamr@2: typedef random::lagged_fibonacci_01<double, 48, 19937, 9842> lagged_fibonacci19937;
williamr@2: typedef random::lagged_fibonacci_01<double, 48, 23209, 13470> lagged_fibonacci23209;
williamr@2: typedef random::lagged_fibonacci_01<double, 48, 44497, 21034> 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