stc/exp/ml_stc_linux_make_v1.0/include/boost/random/lagged_fibonacci.hpp

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