boost/random/linear_congruential.hpp

   1 /* boost random/linear_congruential.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: linear_congruential.hpp 71018 2011-04-05 21:27:52Z steven_watanabe $
  11  *
  12  * Revision history
  13  *  2001-02-18  moved to individual header files
  14  */
  15
  16 #ifndef BOOST_RANDOM_LINEAR_CONGRUENTIAL_HPP
  17 #define BOOST_RANDOM_LINEAR_CONGRUENTIAL_HPP
  18
  19 #include <iostream>
  20 #include <stdexcept>
  21 #include <boost/assert.hpp>
  22 #include <boost/config.hpp>
  23 #include <boost/cstdint.hpp>
  24 #include <boost/limits.hpp>
  25 #include <boost/static_assert.hpp>
  26 #include <boost/integer/static_log2.hpp>
  27 #include <boost/mpl/if.hpp>
  28 #include <boost/type_traits/is_arithmetic.hpp>
  29 #include <boost/random/detail/config.hpp>
  30 #include <boost/random/detail/const_mod.hpp>
  31 #include <boost/random/detail/seed.hpp>
  32 #include <boost/random/detail/seed_impl.hpp>
  33 #include <boost/detail/workaround.hpp>
  34
  35 #include <boost/random/detail/disable_warnings.hpp>
  36
  37 namespace boost {
  38 namespace random {
  39
  40 /**
  41  * Instantiations of class template linear_congruential_engine model a
  42  * \pseudo_random_number_generator. Linear congruential pseudo-random
  43  * number generators are described in:
  44  *
  45  *  @blockquote
  46  *  "Mathematical methods in large-scale computing units", D. H. Lehmer,
  47  *  Proc. 2nd Symposium on Large-Scale Digital Calculating Machines,
  48  *  Harvard University Press, 1951, pp. 141-146
  49  *  @endblockquote
  50  *
  51  * Let x(n) denote the sequence of numbers returned by some pseudo-random
  52  * number generator. Then for the linear congruential generator,
  53  * x(n+1) := (a * x(n) + c) mod m. Parameters for the generator are
  54  * x(0), a, c, m. The template parameter IntType shall denote an integral
  55  * type. It must be large enough to hold values a, c, and m. The template
  56  * parameters a and c must be smaller than m.
  57  *
  58  * Note: The quality of the generator crucially depends on the choice of
  59  * the parameters. User code should use one of the sensibly parameterized
  60  * generators such as minstd_rand instead.
  61  */
  62 template<class IntType, IntType a, IntType c, IntType m>
  63 class linear_congruential_engine
  64 {
  65 public:
  66     typedef IntType result_type;
  67
  68     // Required for old Boost.Random concept
  69     BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
  70
  71     BOOST_STATIC_CONSTANT(IntType, multiplier = a);
  72     BOOST_STATIC_CONSTANT(IntType, increment = c);
  73     BOOST_STATIC_CONSTANT(IntType, modulus = m);
  74     BOOST_STATIC_CONSTANT(IntType, default_seed = 1);
  75
  76     BOOST_STATIC_ASSERT(std::numeric_limits<IntType>::is_integer);
  77     BOOST_STATIC_ASSERT(m == 0 || a < m);
  78     BOOST_STATIC_ASSERT(m == 0 || c < m);
  79
  80     /**
  81      * Constructs a @c linear_congruential_engine, using the default seed
  82      */
  83     linear_congruential_engine() { seed(); }
  84
  85     /**
  86      * Constructs a @c linear_congruential_engine, seeding it with @c x0.
  87      */
  88     BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(linear_congruential_engine,
  89                                                IntType, x0)
  90     { seed(x0); }
  91
  92     /**
  93      * Constructs a @c linear_congruential_engine, seeding it with values
  94      * produced by a call to @c seq.generate().
  95      */
  96     BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(linear_congruential_engine,
  97                                              SeedSeq, seq)
  98     { seed(seq); }
  99
 100     /**
 101      * Constructs a @c linear_congruential_engine  and seeds it
 102      * with values taken from the itrator range [first, last)
 103      * and adjusts first to point to the element after the last one
 104      * used.  If there are not enough elements, throws @c std::invalid_argument.
 105      *
 106      * first and last must be input iterators.
 107      */
 108     template<class It>
 109     linear_congruential_engine(It& first, It last)
 110     {
 111         seed(first, last);
 112     }
 113
 114     // compiler-generated copy constructor and assignment operator are fine
 115
 116     /**
 117      * Calls seed(default_seed)
 118      */
 119     void seed() { seed(default_seed); }
 120
 121     /**
 122      * If c mod m is zero and x0 mod m is zero, changes the current value of
 123      * the generator to 1. Otherwise, changes it to x0 mod m. If c is zero,
 124      * distinct seeds in the range [1,m) will leave the generator in distinct
 125      * states. If c is not zero, the range is [0,m).
 126      */
 127     BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(linear_congruential_engine, IntType, x0)
 128     {
 129         // wrap _x if it doesn't fit in the destination
 130         if(modulus == 0) {
 131             _x = x0;
 132         } else {
 133             _x = x0 % modulus;
 134         }
 135         // handle negative seeds
 136         if(_x <= 0 && _x != 0) {
 137             _x += modulus;
 138         }
 139         // adjust to the correct range
 140         if(increment == 0 && _x == 0) {
 141             _x = 1;
 142         }
 143         BOOST_ASSERT(_x >= (min)());
 144         BOOST_ASSERT(_x <= (max)());
 145     }
 146
 147     /**
 148      * Seeds a @c linear_congruential_engine using values from a SeedSeq.
 149      */
 150     BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(linear_congruential_engine, SeedSeq, seq)
 151     { seed(detail::seed_one_int<IntType, m>(seq)); }
 152
 153     /**
 154      * seeds a @c linear_congruential_engine with values taken
 155      * from the itrator range [first, last) and adjusts @c first to
 156      * point to the element after the last one used.  If there are
 157      * not enough elements, throws @c std::invalid_argument.
 158      *
 159      * @c first and @c last must be input iterators.
 160      */
 161     template<class It>
 162     void seed(It& first, It last)
 163     { seed(detail::get_one_int<IntType, m>(first, last)); }
 164
 165     /**
 166      * Returns the smallest value that the @c linear_congruential_engine
 167      * can produce.
 168      */
 169     static result_type min BOOST_PREVENT_MACRO_SUBSTITUTION ()
 170     { return c == 0 ? 1 : 0; }
 171     /**
 172      * Returns the largest value that the @c linear_congruential_engine
 173      * can produce.
 174      */
 175     static result_type max BOOST_PREVENT_MACRO_SUBSTITUTION ()
 176     { return modulus-1; }
 177
 178     /** Returns the next value of the @c linear_congruential_engine. */
 179     IntType operator()()
 180     {
 181         _x = const_mod<IntType, m>::mult_add(a, _x, c);
 182         return _x;
 183     }
 184
 185     /** Fills a range with random values */
 186     template<class Iter>
 187     void generate(Iter first, Iter last)
 188     { detail::generate_from_int(*this, first, last); }
 189
 190     /** Advances the state of the generator by @c z. */
 191     void discard(boost::uintmax_t z)
 192     {
 193         typedef const_mod<IntType, m> mod_type;
 194         IntType b_inv = mod_type::invert(a-1);
 195         IntType b_gcd = mod_type::mult(a-1, b_inv);
 196         if(b_gcd == 1) {
 197             IntType a_z = mod_type::pow(a, z);
 198             _x = mod_type::mult_add(a_z, _x,
 199                 mod_type::mult(mod_type::mult(c, b_inv), a_z - 1));
 200         } else {
 201             // compute (a^z - 1)*c % (b_gcd * m) / (b / b_gcd) * inv(b / b_gcd)
 202             // we're storing the intermediate result / b_gcd
 203             IntType a_zm1_over_gcd = 0;
 204             IntType a_km1_over_gcd = (a - 1) / b_gcd;
 205             boost::uintmax_t exponent = z;
 206             while(exponent != 0) {
 207                 if(exponent % 2 == 1) {
 208                     a_zm1_over_gcd =
 209                         mod_type::mult_add(
 210                             b_gcd,
 211                             mod_type::mult(a_zm1_over_gcd, a_km1_over_gcd),
 212                             mod_type::add(a_zm1_over_gcd, a_km1_over_gcd));
 213                 }
 214                 a_km1_over_gcd = mod_type::mult_add(
 215                     b_gcd,
 216                     mod_type::mult(a_km1_over_gcd, a_km1_over_gcd),
 217                     mod_type::add(a_km1_over_gcd, a_km1_over_gcd));
 218                 exponent /= 2;
 219             }
 220
 221             IntType a_z = mod_type::mult_add(b_gcd, a_zm1_over_gcd, 1);
 222             IntType num = mod_type::mult(c, a_zm1_over_gcd);
 223             b_inv = mod_type::invert((a-1)/b_gcd);
 224             _x = mod_type::mult_add(a_z, _x, mod_type::mult(b_inv, num));
 225         }
 226     }
 227
 228     friend bool operator==(const linear_congruential_engine& x,
 229                            const linear_congruential_engine& y)
 230     { return x._x == y._x; }
 231     friend bool operator!=(const linear_congruential_engine& x,
 232                            const linear_congruential_engine& y)
 233     { return !(x == y); }
 234
 235 #if !defined(BOOST_RANDOM_NO_STREAM_OPERATORS)
 236     /** Writes a @c linear_congruential_engine to a @c std::ostream. */
 237     template<class CharT, class Traits>
 238     friend std::basic_ostream<CharT,Traits>&
 239     operator<<(std::basic_ostream<CharT,Traits>& os,
 240                const linear_congruential_engine& lcg)
 241     {
 242         return os << lcg._x;
 243     }
 244
 245     /** Reads a @c linear_congruential_engine from a @c std::istream. */
 246     template<class CharT, class Traits>
 247     friend std::basic_istream<CharT,Traits>&
 248     operator>>(std::basic_istream<CharT,Traits>& is,
 249                linear_congruential_engine& lcg)
 250     {
 251         lcg.read(is);
 252         return is;
 253     }
 254 #endif
 255
 256 private:
 257
 258     /// \cond show_private
 259
 260     template<class CharT, class Traits>
 261     void read(std::basic_istream<CharT, Traits>& is) {
 262         IntType x;
 263         if(is >> x) {
 264             if(x >= (min)() && x <= (max)()) {
 265                 _x = x;
 266             } else {
 267                 is.setstate(std::ios_base::failbit);
 268             }
 269         }
 270     }
 271
 272     /// \endcond
 273
 274     IntType _x;
 275 };
 276
 277 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
 278 //  A definition is required even for integral static constants
 279 template<class IntType, IntType a, IntType c, IntType m>
 280 const bool linear_congruential_engine<IntType, a, c, m>::has_fixed_range;
 281 template<class IntType, IntType a, IntType c, IntType m>
 282 const IntType linear_congruential_engine<IntType,a,c,m>::multiplier;
 283 template<class IntType, IntType a, IntType c, IntType m>
 284 const IntType linear_congruential_engine<IntType,a,c,m>::increment;
 285 template<class IntType, IntType a, IntType c, IntType m>
 286 const IntType linear_congruential_engine<IntType,a,c,m>::modulus;
 287 template<class IntType, IntType a, IntType c, IntType m>
 288 const IntType linear_congruential_engine<IntType,a,c,m>::default_seed;
 289 #endif
 290
 291 /// \cond show_deprecated
 292
 293 // provided for backwards compatibility
 294 template<class IntType, IntType a, IntType c, IntType m, IntType val = 0>
 295 class linear_congruential : public linear_congruential_engine<IntType, a, c, m>
 296 {
 297     typedef linear_congruential_engine<IntType, a, c, m> base_type;
 298 public:
 299     linear_congruential(IntType x0 = 1) : base_type(x0) {}
 300     template<class It>
 301     linear_congruential(It& first, It last) : base_type(first, last) {}
 302 };
 303
 304 /// \endcond
 305
 306 /**
 307  * The specialization \minstd_rand0 was originally suggested in
 308  *
 309  *  @blockquote
 310  *  A pseudo-random number generator for the System/360, P.A. Lewis,
 311  *  A.S. Goodman, J.M. Miller, IBM Systems Journal, Vol. 8, No. 2,
 312  *  1969, pp. 136-146
 313  *  @endblockquote
 314  *
 315  * It is examined more closely together with \minstd_rand in
 316  *
 317  *  @blockquote
 318  *  "Random Number Generators: Good ones are hard to find",
 319  *  Stephen K. Park and Keith W. Miller, Communications of
 320  *  the ACM, Vol. 31, No. 10, October 1988, pp. 1192-1201
 321  *  @endblockquote
 322  */
 323 typedef linear_congruential_engine<uint32_t, 16807, 0, 2147483647> minstd_rand0;
 324
 325 /** The specialization \minstd_rand was suggested in
 326  *
 327  *  @blockquote
 328  *  "Random Number Generators: Good ones are hard to find",
 329  *  Stephen K. Park and Keith W. Miller, Communications of
 330  *  the ACM, Vol. 31, No. 10, October 1988, pp. 1192-1201
 331  *  @endblockquote
 332  */
 333 typedef linear_congruential_engine<uint32_t, 48271, 0, 2147483647> minstd_rand;
 334
 335
 336 #if !defined(BOOST_NO_INT64_T) && !defined(BOOST_NO_INTEGRAL_INT64_T)
 337 /**
 338  * Class @c rand48 models a \pseudo_random_number_generator. It uses
 339  * the linear congruential algorithm with the parameters a = 0x5DEECE66D,
 340  * c = 0xB, m = 2**48. It delivers identical results to the @c lrand48()
 341  * function available on some systems (assuming lcong48 has not been called).
 342  *
 343  * It is only available on systems where @c uint64_t is provided as an
 344  * integral type, so that for example static in-class constants and/or
 345  * enum definitions with large @c uint64_t numbers work.
 346  */
 347 class rand48
 348 {
 349 public:
 350     typedef boost::uint32_t result_type;
 351
 352     BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
 353     /**
 354      * Returns the smallest value that the generator can produce
 355      */
 356     static uint32_t min BOOST_PREVENT_MACRO_SUBSTITUTION () { return 0; }
 357     /**
 358      * Returns the largest value that the generator can produce
 359      */
 360     static uint32_t max BOOST_PREVENT_MACRO_SUBSTITUTION ()
 361     { return 0x7FFFFFFF; }
 362
 363     /** Seeds the generator with the default seed. */
 364     rand48() : lcf(cnv(static_cast<uint32_t>(1))) {}
 365     /**
 366      * Constructs a \rand48 generator with x(0) := (x0 << 16) | 0x330e.
 367      */
 368     BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(rand48, result_type, x0)
 369     { seed(x0); }
 370     /**
 371      * Seeds the generator with values produced by @c seq.generate().
 372      */
 373     BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(rand48, SeedSeq, seq)
 374     { seed(seq); }
 375     /**
 376      * Seeds the generator using values from an iterator range,
 377      * and updates first to point one past the last value consumed.
 378      */
 379     template<class It> rand48(It& first, It last) : lcf(first, last) { }
 380
 381     // compiler-generated copy ctor and assignment operator are fine
 382
 383     /** Seeds the generator with the default seed. */
 384     void seed() { seed(static_cast<uint32_t>(1)); }
 385     /**
 386      * Changes the current value x(n) of the generator to (x0 << 16) | 0x330e.
 387      */
 388     BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(rand48, result_type, x0)
 389     { lcf.seed(cnv(x0)); }
 390     /**
 391      * Seeds the generator using values from an iterator range,
 392      * and updates first to point one past the last value consumed.
 393      */
 394     template<class It> void seed(It& first, It last) { lcf.seed(first,last); }
 395     /**
 396      * Seeds the generator with values produced by @c seq.generate().
 397      */
 398     BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(rand48, SeedSeq, seq)
 399     { lcf.seed(seq); }
 400
 401     /**  Returns the next value of the generator. */
 402     uint32_t operator()() { return static_cast<uint32_t>(lcf() >> 17); }
 403
 404     /** Advances the state of the generator by @c z. */
 405     void discard(boost::uintmax_t z) { lcf.discard(z); }
 406
 407     /** Fills a range with random values */
 408     template<class Iter>
 409     void generate(Iter first, Iter last)
 410     {
 411         for(; first != last; ++first) {
 412             *first = (*this)();
 413         }
 414     }
 415
 416 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
 417     /**  Writes a @c rand48 to a @c std::ostream. */
 418     template<class CharT,class Traits>
 419     friend std::basic_ostream<CharT,Traits>&
 420     operator<<(std::basic_ostream<CharT,Traits>& os, const rand48& r)
 421     { os << r.lcf; return os; }
 422
 423     /** Reads a @c rand48 from a @c std::istream. */
 424     template<class CharT,class Traits>
 425     friend std::basic_istream<CharT,Traits>&
 426     operator>>(std::basic_istream<CharT,Traits>& is, rand48& r)
 427     { is >> r.lcf; return is; }
 428 #endif
 429
 430     /**
 431      * Returns true if the two generators will produce identical
 432      * sequences of values.
 433      */
 434     friend bool operator==(const rand48& x, const rand48& y)
 435     { return x.lcf == y.lcf; }
 436     /**
 437      * Returns true if the two generators will produce different
 438      * sequences of values.
 439      */
 440     friend bool operator!=(const rand48& x, const rand48& y)
 441     { return !(x == y); }
 442 private:
 443     /// \cond show_private
 444     typedef random::linear_congruential_engine<uint64_t,
 445         // xxxxULL is not portable
 446         uint64_t(0xDEECE66DUL) | (uint64_t(0x5) << 32),
 447         0xB, uint64_t(1)<<48> lcf_t;
 448     lcf_t lcf;
 449
 450     static boost::uint64_t cnv(boost::uint32_t x)
 451     { return (static_cast<uint64_t>(x) << 16) | 0x330e; }
 452     /// \endcond
 453 };
 454 #endif /* !BOOST_NO_INT64_T && !BOOST_NO_INTEGRAL_INT64_T */
 455
 456 } // namespace random
 457
 458 using random::minstd_rand0;
 459 using random::minstd_rand;
 460 using random::rand48;
 461
 462 } // namespace boost
 463
 464 #include <boost/random/detail/enable_warnings.hpp>
 465
 466 #endif // BOOST_RANDOM_LINEAR_CONGRUENTIAL_HPP