boost/random/subtract_with_carry.hpp

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