boost/random/lagged_fibonacci.hpp

   1 /* boost random/lagged_fibonacci.hpp header file
   2  *
   3  * Copyright Jens Maurer 2000-2001
   4  * Distributed under the Boost Software License, Version 1.0. (See
   5  * accompanying file LICENSE_1_0.txt or copy at
   6  * http://www.boost.org/LICENSE_1_0.txt)
   7  *
   8  * See http://www.boost.org for most recent version including documentation.
   9  *
  10  * $Id$
  11  *
  12  * Revision history
  13  *  2001-02-18  moved to individual header files
  14  */
  15
  16 #ifndef BOOST_RANDOM_LAGGED_FIBONACCI_HPP
  17 #define BOOST_RANDOM_LAGGED_FIBONACCI_HPP
  18
  19 #include <boost/config/no_tr1/cmath.hpp>
  20 #include <iostream>
  21 #include <algorithm>     // std::max
  22 #include <iterator>
  23 #include <boost/config/no_tr1/cmath.hpp>         // std::pow
  24 #include <boost/config.hpp>
  25 #include <boost/limits.hpp>
  26 #include <boost/cstdint.hpp>
  27 #include <boost/detail/workaround.hpp>
  28 #include <boost/random/linear_congruential.hpp>
  29 #include <boost/random/uniform_01.hpp>
  30 #include <boost/random/detail/config.hpp>
  31 #include <boost/random/detail/pass_through_engine.hpp>
  32
  33 namespace boost {
  34 namespace random {
  35
  36 #if BOOST_WORKAROUND(_MSC_FULL_VER, BOOST_TESTED_AT(13102292)) && BOOST_MSVC > 1300
  37 #  define BOOST_RANDOM_EXTRACT_LF
  38 #endif
  39
  40 #if defined(__APPLE_CC__) && defined(__GNUC__) && (__GNUC__ == 3) && (__GNUC_MINOR__ <= 3)
  41 #  define BOOST_RANDOM_EXTRACT_LF
  42 #endif
  43
  44 #  ifdef BOOST_RANDOM_EXTRACT_LF
  45 namespace detail
  46 {
  47   template<class IStream, class F, class RealType>
  48   IStream&
  49   extract_lagged_fibonacci_01(
  50       IStream& is
  51       , F const& f
  52       , unsigned int& i
  53       , RealType* x
  54       , RealType modulus)
  55   {
  56         is >> i >> std::ws;
  57         for(unsigned int i = 0; i < f.long_lag; ++i)
  58         {
  59             RealType value;
  60             is >> value >> std::ws;
  61             x[i] = value / modulus;
  62         }
  63         return is;
  64   }
  65
  66   template<class IStream, class F, class UIntType>
  67   IStream&
  68   extract_lagged_fibonacci(
  69       IStream& is
  70       , F const& f
  71       , unsigned int& i
  72       , UIntType* x)
  73   {
  74       is >> i >> std::ws;
  75       for(unsigned int i = 0; i < f.long_lag; ++i)
  76           is >> x[i] >> std::ws;
  77       return is;
  78   }
  79 }
  80 #  endif
  81
  82 template<class UIntType, int w, unsigned int p, unsigned int q,
  83          UIntType val = 0>
  84 class lagged_fibonacci
  85 {
  86 public:
  87   typedef UIntType result_type;
  88   BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
  89   BOOST_STATIC_CONSTANT(int, word_size = w);
  90   BOOST_STATIC_CONSTANT(unsigned int, long_lag = p);
  91   BOOST_STATIC_CONSTANT(unsigned int, short_lag = q);
  92
  93   result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; }
  94   result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return wordmask; }
  95
  96   lagged_fibonacci() { init_wordmask(); seed(); }
  97   explicit lagged_fibonacci(uint32_t value) { init_wordmask(); seed(value); }
  98   template<class It> lagged_fibonacci(It& first, It last)
  99   { init_wordmask(); seed(first, last); }
 100   // compiler-generated copy ctor and assignment operator are fine
 101
 102 private:
 103   void init_wordmask()
 104   {
 105     wordmask = 0;
 106     for(int j = 0; j < w; ++j)
 107       wordmask |= (1u << j);
 108   }
 109
 110 public:
 111   void seed(uint32_t value = 331u)
 112   {
 113     minstd_rand0 gen(value);
 114     for(unsigned int j = 0; j < long_lag; ++j)
 115       x[j] = gen() & wordmask;
 116     i = long_lag;
 117   }
 118
 119   template<class It>
 120   void seed(It& first, It last)
 121   {
 122     // word size could be smaller than the seed values
 123     unsigned int j;
 124     for(j = 0; j < long_lag && first != last; ++j, ++first)
 125       x[j] = *first & wordmask;
 126     i = long_lag;
 127     if(first == last && j < long_lag)
 128       throw std::invalid_argument("lagged_fibonacci::seed");
 129   }
 130
 131   result_type operator()()
 132   {
 133     if(i >= long_lag)
 134       fill();
 135     return x[i++];
 136   }
 137
 138   static bool validation(result_type x)
 139   {
 140     return x == val;
 141   }
 142
 143 #ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
 144
 145 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
 146   template<class CharT, class Traits>
 147   friend std::basic_ostream<CharT,Traits>&
 148   operator<<(std::basic_ostream<CharT,Traits>& os, const lagged_fibonacci& f)
 149   {
 150     os << f.i << " ";
 151     for(unsigned int i = 0; i < f.long_lag; ++i)
 152       os << f.x[i] << " ";
 153     return os;
 154   }
 155
 156   template<class CharT, class Traits>
 157   friend std::basic_istream<CharT, Traits>&
 158   operator>>(std::basic_istream<CharT, Traits>& is, lagged_fibonacci& f)
 159   {
 160 # ifdef BOOST_RANDOM_EXTRACT_LF
 161       return detail::extract_lagged_fibonacci(is, f, f.i, f.x);
 162 # else
 163       is >> f.i >> std::ws;
 164       for(unsigned int i = 0; i < f.long_lag; ++i)
 165           is >> f.x[i] >> std::ws;
 166       return is;
 167 # endif
 168   }
 169 #endif
 170
 171   friend bool operator==(const lagged_fibonacci& x, const lagged_fibonacci& y)
 172   { return x.i == y.i && std::equal(x.x, x.x+long_lag, y.x); }
 173   friend bool operator!=(const lagged_fibonacci& x,
 174                          const lagged_fibonacci& y)
 175   { return !(x == y); }
 176 #else
 177   // Use a member function; Streamable concept not supported.
 178   bool operator==(const lagged_fibonacci& rhs) const
 179   { return i == rhs.i && std::equal(x, x+long_lag, rhs.x); }
 180   bool operator!=(const lagged_fibonacci& rhs) const
 181   { return !(*this == rhs); }
 182 #endif
 183
 184 private:
 185   void fill();
 186   UIntType wordmask;
 187   unsigned int i;
 188   UIntType x[long_lag];
 189 };
 190
 191 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
 192 //  A definition is required even for integral static constants
 193 template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
 194 const bool lagged_fibonacci<UIntType, w, p, q, val>::has_fixed_range;
 195 template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
 196 const unsigned int lagged_fibonacci<UIntType, w, p, q, val>::long_lag;
 197 template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
 198 const unsigned int lagged_fibonacci<UIntType, w, p, q, val>::short_lag;
 199 #endif
 200
 201 template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
 202 void lagged_fibonacci<UIntType, w, p, q, val>::fill()
 203 {
 204   // two loops to avoid costly modulo operations
 205   {  // extra scope for MSVC brokenness w.r.t. for scope
 206   for(unsigned int j = 0; j < short_lag; ++j)
 207     x[j] = (x[j] + x[j+(long_lag-short_lag)]) & wordmask;
 208   }
 209   for(unsigned int j = short_lag; j < long_lag; ++j)
 210     x[j] = (x[j] + x[j-short_lag]) & wordmask;
 211   i = 0;
 212 }
 213
 214
 215
 216 // lagged Fibonacci generator for the range [0..1)
 217 // contributed by Matthias Troyer
 218 // for p=55, q=24 originally by G. J. Mitchell and D. P. Moore 1958
 219
 220 template<class T, unsigned int p, unsigned int q>
 221 struct fibonacci_validation
 222 {
 223   BOOST_STATIC_CONSTANT(bool, is_specialized = false);
 224   static T value() { return 0; }
 225   static T tolerance() { return 0; }
 226 };
 227
 228 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
 229 //  A definition is required even for integral static constants
 230 template<class T, unsigned int p, unsigned int q>
 231 const bool fibonacci_validation<T, p, q>::is_specialized;
 232 #endif
 233
 234 #define BOOST_RANDOM_FIBONACCI_VAL(T,P,Q,V,E) \
 235 template<> \
 236 struct fibonacci_validation<T, P, Q>  \
 237 {                                     \
 238   BOOST_STATIC_CONSTANT(bool, is_specialized = true);     \
 239   static T value() { return V; }      \
 240   static T tolerance()                \
 241 { return (std::max)(E, static_cast<T>(5*std::numeric_limits<T>::epsilon())); } \
 242 };
 243 // (The extra static_cast<T> in the std::max call above is actually
 244 // unnecessary except for HP aCC 1.30, which claims that
 245 // numeric_limits<double>::epsilon() doesn't actually return a double.)
 246
 247 BOOST_RANDOM_FIBONACCI_VAL(double, 607, 273, 0.4293817707235914, 1e-14)
 248 BOOST_RANDOM_FIBONACCI_VAL(double, 1279, 418, 0.9421630240437659, 1e-14)
 249 BOOST_RANDOM_FIBONACCI_VAL(double, 2281, 1252, 0.1768114046909004, 1e-14)
 250 BOOST_RANDOM_FIBONACCI_VAL(double, 3217, 576, 0.1956232694868209, 1e-14)
 251 BOOST_RANDOM_FIBONACCI_VAL(double, 4423, 2098, 0.9499762202147172, 1e-14)
 252 BOOST_RANDOM_FIBONACCI_VAL(double, 9689, 5502, 0.05737836943695162, 1e-14)
 253 BOOST_RANDOM_FIBONACCI_VAL(double, 19937, 9842, 0.5076528587449834, 1e-14)
 254 BOOST_RANDOM_FIBONACCI_VAL(double, 23209, 13470, 0.5414473810619185, 1e-14)
 255 BOOST_RANDOM_FIBONACCI_VAL(double, 44497,21034, 0.254135073399297, 1e-14)
 256
 257 #undef BOOST_RANDOM_FIBONACCI_VAL
 258
 259 template<class RealType, int w, unsigned int p, unsigned int q>
 260 class lagged_fibonacci_01
 261 {
 262 public:
 263   typedef RealType result_type;
 264   BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
 265   BOOST_STATIC_CONSTANT(int, word_size = w);
 266   BOOST_STATIC_CONSTANT(unsigned int, long_lag = p);
 267   BOOST_STATIC_CONSTANT(unsigned int, short_lag = q);
 268
 269   lagged_fibonacci_01() { init_modulus(); seed(); }
 270   explicit lagged_fibonacci_01(uint32_t value) { init_modulus(); seed(value); }
 271   template<class Generator>
 272   explicit lagged_fibonacci_01(Generator & gen) { init_modulus(); seed(gen); }
 273   template<class It> lagged_fibonacci_01(It& first, It last)
 274   { init_modulus(); seed(first, last); }
 275   // compiler-generated copy ctor and assignment operator are fine
 276
 277 private:
 278   void init_modulus()
 279   {
 280 #ifndef BOOST_NO_STDC_NAMESPACE
 281     // allow for Koenig lookup
 282     using std::pow;
 283 #endif
 284     _modulus = pow(RealType(2), word_size);
 285   }
 286
 287 public:
 288   void seed(uint32_t value = 331u)
 289   {
 290     minstd_rand0 intgen(value);
 291     seed(intgen);
 292   }
 293
 294   // For GCC, moving this function out-of-line prevents inlining, which may
 295   // reduce overall object code size.  However, MSVC does not grok
 296   // out-of-line template member functions.
 297   template<class Generator>
 298   void seed(Generator & gen)
 299   {
 300     // use pass-by-reference, but wrap argument in pass_through_engine
 301     typedef detail::pass_through_engine<Generator&> ref_gen;
 302     uniform_01<ref_gen, RealType> gen01 =
 303       uniform_01<ref_gen, RealType>(ref_gen(gen));
 304     // I could have used std::generate_n, but it takes "gen" by value
 305     for(unsigned int j = 0; j < long_lag; ++j)
 306       x[j] = gen01();
 307     i = long_lag;
 308   }
 309
 310   template<class It>
 311   void seed(It& first, It last)
 312   {
 313 #ifndef BOOST_NO_STDC_NAMESPACE
 314     // allow for Koenig lookup
 315     using std::fmod;
 316     using std::pow;
 317 #endif
 318     unsigned long mask = ~((~0u) << (w%32));   // now lowest w bits set
 319     RealType two32 = pow(RealType(2), 32);
 320     unsigned int j;
 321     for(j = 0; j < long_lag && first != last; ++j) {
 322       x[j] = RealType(0);
 323       for(int k = 0; k < w/32 && first != last; ++k, ++first)
 324         x[j] += *first / pow(two32,k+1);
 325       if(first != last && mask != 0) {
 326         x[j] += fmod((*first & mask) / _modulus, RealType(1));
 327         ++first;
 328       }
 329     }
 330     i = long_lag;
 331     if(first == last && j < long_lag)
 332       throw std::invalid_argument("lagged_fibonacci_01::seed");
 333   }
 334
 335   result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(0); }
 336   result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(1); }
 337
 338   result_type operator()()
 339   {
 340     if(i >= long_lag)
 341       fill();
 342     return x[i++];
 343   }
 344
 345   static bool validation(result_type x)
 346   {
 347     result_type v = fibonacci_validation<result_type, p, q>::value();
 348     result_type epsilon = fibonacci_validation<result_type, p, q>::tolerance();
 349     // std::abs is a source of trouble: sometimes, it's not overloaded
 350     // for double, plus the usual namespace std noncompliance -> avoid it
 351     // using std::abs;
 352     // return abs(x - v) < 5 * epsilon
 353     return x > v - epsilon && x < v + epsilon;
 354   }
 355
 356 #ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
 357
 358 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
 359   template<class CharT, class Traits>
 360   friend std::basic_ostream<CharT,Traits>&
 361   operator<<(std::basic_ostream<CharT,Traits>& os, const lagged_fibonacci_01&f)
 362   {
 363 #ifndef BOOST_NO_STDC_NAMESPACE
 364     // allow for Koenig lookup
 365     using std::pow;
 366 #endif
 367     os << f.i << " ";
 368     std::ios_base::fmtflags oldflags = os.flags(os.dec | os.fixed | os.left);
 369     for(unsigned int i = 0; i < f.long_lag; ++i)
 370       os << f.x[i] * f._modulus << " ";
 371     os.flags(oldflags);
 372     return os;
 373   }
 374
 375   template<class CharT, class Traits>
 376   friend std::basic_istream<CharT, Traits>&
 377   operator>>(std::basic_istream<CharT, Traits>& is, lagged_fibonacci_01& f)
 378     {
 379 # ifdef BOOST_RANDOM_EXTRACT_LF
 380         return detail::extract_lagged_fibonacci_01(is, f, f.i, f.x, f._modulus);
 381 # else
 382         is >> f.i >> std::ws;
 383         for(unsigned int i = 0; i < f.long_lag; ++i) {
 384             typename lagged_fibonacci_01::result_type value;
 385             is >> value >> std::ws;
 386             f.x[i] = value / f._modulus;
 387         }
 388         return is;
 389 # endif
 390     }
 391 #endif
 392
 393   friend bool operator==(const lagged_fibonacci_01& x,
 394                          const lagged_fibonacci_01& y)
 395   { return x.i == y.i && std::equal(x.x, x.x+long_lag, y.x); }
 396   friend bool operator!=(const lagged_fibonacci_01& x,
 397                          const lagged_fibonacci_01& y)
 398   { return !(x == y); }
 399 #else
 400   // Use a member function; Streamable concept not supported.
 401   bool operator==(const lagged_fibonacci_01& rhs) const
 402   { return i == rhs.i && std::equal(x, x+long_lag, rhs.x); }
 403   bool operator!=(const lagged_fibonacci_01& rhs) const
 404   { return !(*this == rhs); }
 405 #endif
 406
 407 private:
 408   void fill();
 409   unsigned int i;
 410   RealType x[long_lag];
 411   RealType _modulus;
 412 };
 413
 414 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
 415 //  A definition is required even for integral static constants
 416 template<class RealType, int w, unsigned int p, unsigned int q>
 417 const bool lagged_fibonacci_01<RealType, w, p, q>::has_fixed_range;
 418 template<class RealType, int w, unsigned int p, unsigned int q>
 419 const unsigned int lagged_fibonacci_01<RealType, w, p, q>::long_lag;
 420 template<class RealType, int w, unsigned int p, unsigned int q>
 421 const unsigned int lagged_fibonacci_01<RealType, w, p, q>::short_lag;
 422 template<class RealType, int w, unsigned int p, unsigned int q>
 423 const int lagged_fibonacci_01<RealType,w,p,q>::word_size;
 424
 425 #endif
 426
 427 template<class RealType, int w, unsigned int p, unsigned int q>
 428 void lagged_fibonacci_01<RealType, w, p, q>::fill()
 429 {
 430   // two loops to avoid costly modulo operations
 431   {  // extra scope for MSVC brokenness w.r.t. for scope
 432   for(unsigned int j = 0; j < short_lag; ++j) {
 433     RealType t = x[j] + x[j+(long_lag-short_lag)];
 434     if(t >= RealType(1))
 435       t -= RealType(1);
 436     x[j] = t;
 437   }
 438   }
 439   for(unsigned int j = short_lag; j < long_lag; ++j) {
 440     RealType t = x[j] + x[j-short_lag];
 441     if(t >= RealType(1))
 442       t -= RealType(1);
 443     x[j] = t;
 444   }
 445   i = 0;
 446 }
 447
 448 } // namespace random
 449
 450 typedef random::lagged_fibonacci_01<double, 48, 607, 273> lagged_fibonacci607;
 451 typedef random::lagged_fibonacci_01<double, 48, 1279, 418> lagged_fibonacci1279;
 452 typedef random::lagged_fibonacci_01<double, 48, 2281, 1252> lagged_fibonacci2281;
 453 typedef random::lagged_fibonacci_01<double, 48, 3217, 576> lagged_fibonacci3217;
 454 typedef random::lagged_fibonacci_01<double, 48, 4423, 2098> lagged_fibonacci4423;
 455 typedef random::lagged_fibonacci_01<double, 48, 9689, 5502> lagged_fibonacci9689;
 456 typedef random::lagged_fibonacci_01<double, 48, 19937, 9842> lagged_fibonacci19937;
 457 typedef random::lagged_fibonacci_01<double, 48, 23209, 13470> lagged_fibonacci23209;
 458 typedef random::lagged_fibonacci_01<double, 48, 44497, 21034> lagged_fibonacci44497;
 459
 460
 461 // It is possible to partially specialize uniform_01<> on lagged_fibonacci_01<>
 462 // to help the compiler generate efficient code.  For GCC, this seems useless,
 463 // because GCC optimizes (x-0)/(1-0) to (x-0).  This is good enough for now.
 464
 465 } // namespace boost
 466
 467 #endif // BOOST_RANDOM_LAGGED_FIBONACCI_HPP