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