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

   1 /* boost random/subtract_with_carry.hpp header file\r
   2  *\r
   3  * Copyright Jens Maurer 2002\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: subtract_with_carry.hpp 53871 2009-06-13 17:54:06Z steven_watanabe $\r
  11  *\r
  12  * Revision history\r
  13  *  2002-03-02  created\r
  14  */\r
  15 \r
  16 #ifndef BOOST_RANDOM_SUBTRACT_WITH_CARRY_HPP\r
  17 #define BOOST_RANDOM_SUBTRACT_WITH_CARRY_HPP\r
  18 \r
  19 #include <boost/config/no_tr1/cmath.hpp>\r
  20 #include <iostream>\r
  21 #include <algorithm>     // std::equal\r
  22 #include <stdexcept>\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/static_assert.hpp>\r
  28 #include <boost/detail/workaround.hpp>\r
  29 #include <boost/random/detail/config.hpp>\r
  30 #include <boost/random/detail/seed.hpp>\r
  31 #include <boost/random/linear_congruential.hpp>\r
  32 \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_SWC_01\r
  39 #endif\r
  40 \r
  41 #if defined(__APPLE_CC__) && defined(__GNUC__) && (__GNUC__ == 3) && (__GNUC_MINOR__ <= 3)\r
  42 #  define BOOST_RANDOM_EXTRACT_SWC_01\r
  43 #endif\r
  44 \r
  45 # ifdef BOOST_RANDOM_EXTRACT_SWC_01\r
  46 namespace detail\r
  47 {\r
  48   template <class IStream, class SubtractWithCarry, class RealType>\r
  49   void extract_subtract_with_carry_01(\r
  50       IStream& is\r
  51       , SubtractWithCarry& f\r
  52       , RealType& carry\r
  53       , RealType* x\r
  54       , RealType modulus)\r
  55   {\r
  56     RealType value;\r
  57     for(unsigned int j = 0; j < f.long_lag; ++j) {\r
  58       is >> value >> std::ws;\r
  59       x[j] = value / modulus;\r
  60     }\r
  61     is >> value >> std::ws;\r
  62     carry = value / modulus;\r
  63   }\r
  64 }\r
  65 # endif \r
  66 // subtract-with-carry generator\r
  67 // Marsaglia and Zaman\r
  68 \r
  69 template<class IntType, IntType m, unsigned int s, unsigned int r,\r
  70   IntType val>\r
  71 class subtract_with_carry\r
  72 {\r
  73 public:\r
  74   typedef IntType result_type;\r
  75   BOOST_STATIC_CONSTANT(bool, has_fixed_range = true);\r
  76   BOOST_STATIC_CONSTANT(result_type, min_value = 0);\r
  77   BOOST_STATIC_CONSTANT(result_type, max_value = m-1);\r
  78   BOOST_STATIC_CONSTANT(result_type, modulus = m);\r
  79   BOOST_STATIC_CONSTANT(unsigned int, long_lag = r);\r
  80   BOOST_STATIC_CONSTANT(unsigned int, short_lag = s);\r
  81 \r
  82   subtract_with_carry() {\r
  83     // MSVC fails BOOST_STATIC_ASSERT with std::numeric_limits at class scope\r
  84 #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS\r
  85     BOOST_STATIC_ASSERT(std::numeric_limits<result_type>::is_signed);\r
  86     BOOST_STATIC_ASSERT(std::numeric_limits<result_type>::is_integer);\r
  87 #endif\r
  88     seed();\r
  89   }\r
  90   BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(subtract_with_carry, uint32_t, value)\r
  91   { seed(value); }\r
  92   BOOST_RANDOM_DETAIL_GENERATOR_CONSTRUCTOR(subtract_with_carry, Generator, gen)\r
  93   { seed(gen); }\r
  94   template<class It> subtract_with_carry(It& first, It last) { seed(first,last); }\r
  95 \r
  96   // compiler-generated copy ctor and assignment operator are fine\r
  97 \r
  98   void seed() { seed(19780503u); }\r
  99   BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(subtract_with_carry, uint32_t, value)\r
 100   {\r
 101     random::linear_congruential<int32_t, 40014, 0, 2147483563, 0> intgen(value);\r
 102     seed(intgen);\r
 103   }\r
 104 \r
 105   // For GCC, moving this function out-of-line prevents inlining, which may\r
 106   // reduce overall object code size.  However, MSVC does not grok\r
 107   // out-of-line template member functions.\r
 108   BOOST_RANDOM_DETAIL_GENERATOR_SEED(subtract_with_carry, Generator, gen)\r
 109   {\r
 110     // I could have used std::generate_n, but it takes "gen" by value\r
 111     for(unsigned int j = 0; j < long_lag; ++j)\r
 112       x[j] = gen() % modulus;\r
 113     carry = (x[long_lag-1] == 0);\r
 114     k = 0;\r
 115   }\r
 116 \r
 117   template<class It>\r
 118   void seed(It& first, It last)\r
 119   {\r
 120     unsigned int j;\r
 121     for(j = 0; j < long_lag && first != last; ++j, ++first)\r
 122       x[j] = *first % modulus;\r
 123     if(first == last && j < long_lag)\r
 124       throw std::invalid_argument("subtract_with_carry::seed");\r
 125     carry = (x[long_lag-1] == 0);\r
 126     k = 0;\r
 127    }\r
 128 \r
 129   result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return min_value; }\r
 130   result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return max_value; }\r
 131 \r
 132   result_type operator()()\r
 133   {\r
 134     int short_index = k - short_lag;\r
 135     if(short_index < 0)\r
 136       short_index += long_lag;\r
 137     IntType delta;\r
 138     if (x[short_index] >= x[k] + carry) {\r
 139       // x(n) >= 0\r
 140       delta =  x[short_index] - (x[k] + carry);\r
 141       carry = 0;\r
 142     } else {\r
 143       // x(n) < 0\r
 144       delta = modulus - x[k] - carry + x[short_index];\r
 145       carry = 1;\r
 146     }\r
 147     x[k] = delta;\r
 148     ++k;\r
 149     if(k >= long_lag)\r
 150       k = 0;\r
 151     return delta;\r
 152   }\r
 153 \r
 154 public:\r
 155   static bool validation(result_type x) { return x == val; }\r
 156   \r
 157 #ifndef BOOST_NO_OPERATORS_IN_NAMESPACE\r
 158 \r
 159 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS\r
 160   template<class CharT, class Traits>\r
 161   friend std::basic_ostream<CharT,Traits>&\r
 162   operator<<(std::basic_ostream<CharT,Traits>& os,\r
 163              const subtract_with_carry& f)\r
 164   {\r
 165     for(unsigned int j = 0; j < f.long_lag; ++j)\r
 166       os << f.compute(j) << " ";\r
 167     os << f.carry << " ";\r
 168     return os;\r
 169   }\r
 170 \r
 171   template<class CharT, class Traits>\r
 172   friend std::basic_istream<CharT,Traits>&\r
 173   operator>>(std::basic_istream<CharT,Traits>& is, subtract_with_carry& f)\r
 174   {\r
 175     for(unsigned int j = 0; j < f.long_lag; ++j)\r
 176       is >> f.x[j] >> std::ws;\r
 177     is >> f.carry >> std::ws;\r
 178     f.k = 0;\r
 179     return is;\r
 180   }\r
 181 #endif\r
 182 \r
 183   friend bool operator==(const subtract_with_carry& x, const subtract_with_carry& y)\r
 184   {\r
 185     for(unsigned int j = 0; j < r; ++j)\r
 186       if(x.compute(j) != y.compute(j))\r
 187         return false;\r
 188     return true;\r
 189   }\r
 190 \r
 191   friend bool operator!=(const subtract_with_carry& x, const subtract_with_carry& y)\r
 192   { return !(x == y); }\r
 193 #else\r
 194   // Use a member function; Streamable concept not supported.\r
 195   bool operator==(const subtract_with_carry& rhs) const\r
 196   {\r
 197     for(unsigned int j = 0; j < r; ++j)\r
 198       if(compute(j) != rhs.compute(j))\r
 199         return false;\r
 200     return true;\r
 201   }\r
 202 \r
 203   bool operator!=(const subtract_with_carry& rhs) const\r
 204   { return !(*this == rhs); }\r
 205 #endif\r
 206 \r
 207 private:\r
 208   // returns x(i-r+index), where index is in 0..r-1\r
 209   IntType compute(unsigned int index) const\r
 210   {\r
 211     return x[(k+index) % long_lag];\r
 212   }\r
 213 \r
 214   // state representation; next output (state) is x(i)\r
 215   //   x[0]  ... x[k] x[k+1] ... x[long_lag-1]     represents\r
 216   //  x(i-k) ... x(i) x(i+1) ... x(i-k+long_lag-1)\r
 217   // speed: base: 20-25 nsec\r
 218   // ranlux_4: 230 nsec, ranlux_7: 430 nsec, ranlux_14: 810 nsec\r
 219   // This state representation makes operator== and save/restore more\r
 220   // difficult, because we've already computed "too much" and thus\r
 221   // have to undo some steps to get at x(i-r) etc.\r
 222 \r
 223   // state representation: next output (state) is x(i)\r
 224   //   x[0]  ... x[k] x[k+1]          ... x[long_lag-1]     represents\r
 225   //  x(i-k) ... x(i) x(i-long_lag+1) ... x(i-k-1)\r
 226   // speed: base 28 nsec\r
 227   // ranlux_4: 370 nsec, ranlux_7: 688 nsec, ranlux_14: 1343 nsec\r
 228   IntType x[long_lag];\r
 229   unsigned int k;\r
 230   int carry;\r
 231 };\r
 232 \r
 233 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION\r
 234 //  A definition is required even for integral static constants\r
 235 template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val>\r
 236 const bool subtract_with_carry<IntType, m, s, r, val>::has_fixed_range;\r
 237 template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val>\r
 238 const IntType subtract_with_carry<IntType, m, s, r, val>::min_value;\r
 239 template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val>\r
 240 const IntType subtract_with_carry<IntType, m, s, r, val>::max_value;\r
 241 template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val>\r
 242 const IntType subtract_with_carry<IntType, m, s, r, val>::modulus;\r
 243 template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val>\r
 244 const unsigned int subtract_with_carry<IntType, m, s, r, val>::long_lag;\r
 245 template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val>\r
 246 const unsigned int subtract_with_carry<IntType, m, s, r, val>::short_lag;\r
 247 #endif\r
 248 \r
 249 \r
 250 // use a floating-point representation to produce values in [0..1)\r
 251 template<class RealType, int w, unsigned int s, unsigned int r, int val=0>\r
 252 class subtract_with_carry_01\r
 253 {\r
 254 public:\r
 255   typedef RealType result_type;\r
 256   BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);\r
 257   BOOST_STATIC_CONSTANT(int, word_size = w);\r
 258   BOOST_STATIC_CONSTANT(unsigned int, long_lag = r);\r
 259   BOOST_STATIC_CONSTANT(unsigned int, short_lag = s);\r
 260 \r
 261 #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS\r
 262   BOOST_STATIC_ASSERT(!std::numeric_limits<result_type>::is_integer);\r
 263 #endif\r
 264 \r
 265   subtract_with_carry_01() { init_modulus(); seed(); }\r
 266   explicit subtract_with_carry_01(uint32_t value)\r
 267   { init_modulus(); seed(value);   }\r
 268   template<class It> subtract_with_carry_01(It& first, It last)\r
 269   { init_modulus(); seed(first,last); }\r
 270 \r
 271 private:\r
 272   void init_modulus()\r
 273   {\r
 274 #ifndef BOOST_NO_STDC_NAMESPACE\r
 275     // allow for Koenig lookup\r
 276     using std::pow;\r
 277 #endif\r
 278     _modulus = pow(RealType(2), word_size);\r
 279   }\r
 280 \r
 281 public:\r
 282   // compiler-generated copy ctor and assignment operator are fine\r
 283 \r
 284   void seed(uint32_t value = 19780503u)\r
 285   {\r
 286 #ifndef BOOST_NO_STDC_NAMESPACE\r
 287     // allow for Koenig lookup\r
 288     using std::fmod;\r
 289 #endif\r
 290     random::linear_congruential<int32_t, 40014, 0, 2147483563, 0> gen(value);\r
 291     unsigned long array[(w+31)/32 * long_lag];\r
 292     for(unsigned int j = 0; j < sizeof(array)/sizeof(unsigned long); ++j)\r
 293       array[j] = gen();\r
 294     unsigned long * start = array;\r
 295     seed(start, array + sizeof(array)/sizeof(unsigned long));\r
 296   }\r
 297 \r
 298   template<class It>\r
 299   void seed(It& first, It last)\r
 300   {\r
 301 #ifndef BOOST_NO_STDC_NAMESPACE\r
 302     // allow for Koenig lookup\r
 303     using std::fmod;\r
 304     using std::pow;\r
 305 #endif\r
 306     unsigned long mask = ~((~0u) << (w%32));   // now lowest (w%32) bits set\r
 307     RealType two32 = pow(RealType(2), 32);\r
 308     unsigned int j;\r
 309     for(j = 0; j < long_lag && first != last; ++j) {\r
 310       x[j] = RealType(0);\r
 311       for(int i = 0; i < w/32 && first != last; ++i, ++first)\r
 312         x[j] += *first / pow(two32,i+1);\r
 313       if(first != last && mask != 0) {\r
 314         x[j] += fmod((*first & mask) / _modulus, RealType(1));\r
 315         ++first;\r
 316       }\r
 317     }\r
 318     if(first == last && j < long_lag)\r
 319       throw std::invalid_argument("subtract_with_carry_01::seed");\r
 320     carry = (x[long_lag-1] ? 0 : 1 / _modulus);\r
 321     k = 0;\r
 322   }\r
 323 \r
 324   result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(0); }\r
 325   result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(1); }\r
 326 \r
 327   result_type operator()()\r
 328   {\r
 329     int short_index = k - short_lag;\r
 330     if(short_index < 0)\r
 331       short_index += long_lag;\r
 332     RealType delta = x[short_index] - x[k] - carry;\r
 333     if(delta < 0) {\r
 334       delta += RealType(1);\r
 335       carry = RealType(1)/_modulus;\r
 336     } else {\r
 337       carry = 0;\r
 338     }\r
 339     x[k] = delta;\r
 340     ++k;\r
 341     if(k >= long_lag)\r
 342       k = 0;\r
 343     return delta;\r
 344   }\r
 345 \r
 346   static bool validation(result_type x)\r
 347   { return x == val/pow(RealType(2), word_size); }\r
 348   \r
 349 #ifndef BOOST_NO_OPERATORS_IN_NAMESPACE\r
 350 \r
 351 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS\r
 352   template<class CharT, class Traits>\r
 353   friend std::basic_ostream<CharT,Traits>&\r
 354   operator<<(std::basic_ostream<CharT,Traits>& os,\r
 355              const subtract_with_carry_01& f)\r
 356   {\r
 357 #ifndef BOOST_NO_STDC_NAMESPACE\r
 358     // allow for Koenig lookup\r
 359     using std::pow;\r
 360 #endif\r
 361     std::ios_base::fmtflags oldflags = os.flags(os.dec | os.fixed | os.left); \r
 362     for(unsigned int j = 0; j < f.long_lag; ++j)\r
 363       os << (f.compute(j) * f._modulus) << " ";\r
 364     os << (f.carry * f._modulus);\r
 365     os.flags(oldflags);\r
 366     return os;\r
 367   }\r
 368 \r
 369   template<class CharT, class Traits>\r
 370   friend std::basic_istream<CharT,Traits>&\r
 371   operator>>(std::basic_istream<CharT,Traits>& is, subtract_with_carry_01& f)\r
 372   {\r
 373 # ifdef BOOST_RANDOM_EXTRACT_SWC_01\r
 374       detail::extract_subtract_with_carry_01(is, f, f.carry, f.x, f._modulus);\r
 375 # else\r
 376     // MSVC (up to 7.1) and Borland (up to 5.64) don't handle the template type\r
 377     // parameter "RealType" available from the class template scope, so use\r
 378     // the member typedef\r
 379     typename subtract_with_carry_01::result_type value;\r
 380     for(unsigned int j = 0; j < long_lag; ++j) {\r
 381       is >> value >> std::ws;\r
 382       f.x[j] = value / f._modulus;\r
 383     }\r
 384     is >> value >> std::ws;\r
 385     f.carry = value / f._modulus;\r
 386 # endif \r
 387     f.k = 0;\r
 388     return is;\r
 389   }\r
 390 #endif\r
 391 \r
 392   friend bool operator==(const subtract_with_carry_01& x,\r
 393                          const subtract_with_carry_01& y)\r
 394   {\r
 395     for(unsigned int j = 0; j < r; ++j)\r
 396       if(x.compute(j) != y.compute(j))\r
 397         return false;\r
 398     return true;\r
 399   }\r
 400 \r
 401   friend bool operator!=(const subtract_with_carry_01& x,\r
 402                          const subtract_with_carry_01& y)\r
 403   { return !(x == y); }\r
 404 #else\r
 405   // Use a member function; Streamable concept not supported.\r
 406   bool operator==(const subtract_with_carry_01& rhs) const\r
 407   { \r
 408     for(unsigned int j = 0; j < r; ++j)\r
 409       if(compute(j) != rhs.compute(j))\r
 410         return false;\r
 411     return true;\r
 412   }\r
 413 \r
 414   bool operator!=(const subtract_with_carry_01& rhs) const\r
 415   { return !(*this == rhs); }\r
 416 #endif\r
 417 \r
 418 private:\r
 419   RealType compute(unsigned int index) const;\r
 420   unsigned int k;\r
 421   RealType carry;\r
 422   RealType x[long_lag];\r
 423   RealType _modulus;\r
 424 };\r
 425 \r
 426 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION\r
 427 //  A definition is required even for integral static constants\r
 428 template<class RealType, int w, unsigned int s, unsigned int r, int val>\r
 429 const bool subtract_with_carry_01<RealType, w, s, r, val>::has_fixed_range;\r
 430 template<class RealType, int w, unsigned int s, unsigned int r, int val>\r
 431 const int subtract_with_carry_01<RealType, w, s, r, val>::word_size;\r
 432 template<class RealType, int w, unsigned int s, unsigned int r, int val>\r
 433 const unsigned int subtract_with_carry_01<RealType, w, s, r, val>::long_lag;\r
 434 template<class RealType, int w, unsigned int s, unsigned int r, int val>\r
 435 const unsigned int subtract_with_carry_01<RealType, w, s, r, val>::short_lag;\r
 436 #endif\r
 437 \r
 438 template<class RealType, int w, unsigned int s, unsigned int r, int val>\r
 439 RealType subtract_with_carry_01<RealType, w, s, r, val>::compute(unsigned int index) const\r
 440 {\r
 441   return x[(k+index) % long_lag];\r
 442 }\r
 443 \r
 444 \r
 445 } // namespace random\r
 446 } // namespace boost\r
 447 \r
 448 #endif // BOOST_RANDOM_SUBTRACT_WITH_CARRY_HPP\r