boost/random/lagged_fibonacci.hpp

   1 /* boost random/lagged_fibonacci.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: lagged_fibonacci.hpp 72951 2011-07-07 04:57:37Z steven_watanabe $
  11  *
  12  * Revision history
  13  *  2001-02-18  moved to individual header files
  14  */
  15
  16 #ifndef BOOST_RANDOM_LAGGED_FIBONACCI_HPP
  17 #define BOOST_RANDOM_LAGGED_FIBONACCI_HPP
  18
  19 #include <istream>
  20 #include <iosfwd>
  21 #include <algorithm>     // std::max
  22 #include <iterator>
  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/integer/integer_mask.hpp>
  28 #include <boost/random/linear_congruential.hpp>
  29 #include <boost/random/uniform_01.hpp>
  30 #include <boost/random/detail/config.hpp>
  31 #include <boost/random/detail/seed.hpp>
  32 #include <boost/random/detail/operators.hpp>
  33 #include <boost/random/detail/generator_seed_seq.hpp>
  34
  35 namespace boost {
  36 namespace random {
  37
  38 /**
  39  * Instantiations of class template \lagged_fibonacci_engine model a
  40  * \pseudo_random_number_generator. It uses a lagged Fibonacci
  41  * algorithm with two lags @c p and @c q:
  42  * x(i) = x(i-p) + x(i-q) (mod 2<sup>w</sup>) with p > q.
  43  */
  44 template<class UIntType, int w, unsigned int p, unsigned int q>
  45 class lagged_fibonacci_engine
  46 {
  47 public:
  48     typedef UIntType result_type;
  49     BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
  50     BOOST_STATIC_CONSTANT(int, word_size = w);
  51     BOOST_STATIC_CONSTANT(unsigned int, long_lag = p);
  52     BOOST_STATIC_CONSTANT(unsigned int, short_lag = q);
  53
  54     BOOST_STATIC_CONSTANT(UIntType, default_seed = 331u);
  55
  56     /** Returns the smallest value that the generator can produce. */
  57     static result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () { return 0; }
  58     /** Returns the largest value that the generator can produce. */
  59     static result_type max BOOST_PREVENT_MACRO_SUBSTITUTION ()
  60     { return low_bits_mask_t<w>::sig_bits; }
  61
  62     /** Creates a new @c lagged_fibonacci_engine and calls @c seed(). */
  63     lagged_fibonacci_engine() { seed(); }
  64
  65     /** Creates a new @c lagged_fibonacci_engine and calls @c seed(value). */
  66     BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(lagged_fibonacci_engine,
  67         UIntType, value)
  68     { seed(value); }
  69
  70     /** Creates a new @c lagged_fibonacci_engine and calls @c seed(seq). */
  71     BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(lagged_fibonacci_engine,
  72         SeedSeq, seq)
  73     { seed(seq); }
  74
  75     /**
  76      * Creates a new @c lagged_fibonacci_engine and calls @c seed(first, last).
  77      */
  78     template<class It> lagged_fibonacci_engine(It& first, It last)
  79     { seed(first, last); }
  80
  81     // compiler-generated copy ctor and assignment operator are fine
  82
  83     /** Calls @c seed(default_seed). */
  84     void seed() { seed(default_seed); }
  85
  86     /**
  87      * Sets the state of the generator to values produced by
  88      * a \minstd_rand0 generator.
  89      */
  90     BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(lagged_fibonacci_engine,
  91         UIntType, value)
  92     {
  93         minstd_rand0 intgen(static_cast<boost::uint32_t>(value));
  94         detail::generator_seed_seq<minstd_rand0> gen(intgen);
  95         seed(gen);
  96     }
  97
  98     /**
  99      * Sets the state of the generator using values produced by seq.
 100      */
 101     BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(lagged_fibonacci_engine, SeedSeq, seq)
 102     {
 103         detail::seed_array_int<w>(seq, x);
 104         i = long_lag;
 105     }
 106
 107     /**
 108      * Sets the state of the generator to values from the iterator
 109      * range [first, last).  If there are not enough elements in the
 110      * range [first, last) throws @c std::invalid_argument.
 111      */
 112     template<class It>
 113     void seed(It& first, It last)
 114     {
 115         detail::fill_array_int<w>(first, last, x);
 116         i = long_lag;
 117     }
 118
 119     /** Returns the next value of the generator. */
 120     result_type operator()()
 121     {
 122         if(i >= long_lag)
 123             fill();
 124         return x[i++];
 125     }
 126
 127     /** Fills a range with random values */
 128     template<class Iter>
 129     void generate(Iter first, Iter last)
 130     { detail::generate_from_int(*this, first, last); }
 131
 132     /** Advances the state of the generator by @c z. */
 133     void discard(boost::uintmax_t z)
 134     {
 135         for(boost::uintmax_t j = 0; j < z; ++j) {
 136             (*this)();
 137         }
 138     }
 139
 140     /**
 141      * Writes the textual representation of the generator to a @c std::ostream.
 142      */
 143     BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, lagged_fibonacci_engine, f)
 144     {
 145         os << f.i;
 146         for(unsigned int i = 0; i < f.long_lag; ++i)
 147             os << ' ' << f.x[i];
 148         return os;
 149     }
 150
 151     /**
 152      * Reads the textual representation of the generator from a @c std::istream.
 153      */
 154     BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, lagged_fibonacci_engine, f)
 155     {
 156         is >> f.i >> std::ws;
 157         for(unsigned int i = 0; i < f.long_lag; ++i)
 158             is >> f.x[i] >> std::ws;
 159         return is;
 160     }
 161
 162     /**
 163      * Returns true if the two generators will produce identical
 164      * sequences of outputs.
 165      */
 166     BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(lagged_fibonacci_engine, x, y)
 167     { return x.i == y.i && std::equal(x.x, x.x+long_lag, y.x); }
 168
 169     /**
 170      * Returns true if the two generators will produce different
 171      * sequences of outputs.
 172      */
 173     BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(lagged_fibonacci_engine)
 174
 175 private:
 176     /// \cond show_private
 177     void fill();
 178     /// \endcond
 179
 180     unsigned int i;
 181     UIntType x[long_lag];
 182 };
 183
 184 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
 185 //  A definition is required even for integral static constants
 186 template<class UIntType, int w, unsigned int p, unsigned int q>
 187 const bool lagged_fibonacci_engine<UIntType, w, p, q>::has_fixed_range;
 188 template<class UIntType, int w, unsigned int p, unsigned int q>
 189 const unsigned int lagged_fibonacci_engine<UIntType, w, p, q>::long_lag;
 190 template<class UIntType, int w, unsigned int p, unsigned int q>
 191 const unsigned int lagged_fibonacci_engine<UIntType, w, p, q>::short_lag;
 192 template<class UIntType, int w, unsigned int p, unsigned int q>
 193 const UIntType lagged_fibonacci_engine<UIntType, w, p, q>::default_seed;
 194 #endif
 195
 196 /// \cond show_private
 197
 198 template<class UIntType, int w, unsigned int p, unsigned int q>
 199 void lagged_fibonacci_engine<UIntType, w, p, q>::fill()
 200 {
 201     // two loops to avoid costly modulo operations
 202     {  // extra scope for MSVC brokenness w.r.t. for scope
 203     for(unsigned int j = 0; j < short_lag; ++j)
 204         x[j] = (x[j] + x[j+(long_lag-short_lag)]) & low_bits_mask_t<w>::sig_bits;
 205     }
 206     for(unsigned int j = short_lag; j < long_lag; ++j)
 207         x[j] = (x[j] + x[j-short_lag]) & low_bits_mask_t<w>::sig_bits;
 208     i = 0;
 209 }
 210
 211 /// \endcond
 212
 213 /// \cond show_deprecated
 214
 215 // provided for backwards compatibility
 216 template<class UIntType, int w, unsigned int p, unsigned int q, UIntType v = 0>
 217 class lagged_fibonacci : public lagged_fibonacci_engine<UIntType, w, p, q>
 218 {
 219     typedef lagged_fibonacci_engine<UIntType, w, p, q> base_type;
 220 public:
 221     lagged_fibonacci() {}
 222     BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(lagged_fibonacci, UIntType, val)
 223     { this->seed(val); }
 224     BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(lagged_fibonacci, SeedSeq, seq)
 225     { this->seed(seq); }
 226     template<class It>
 227     lagged_fibonacci(It& first, It last) : base_type(first, last) {}
 228 };
 229
 230 /// \endcond
 231
 232 // lagged Fibonacci generator for the range [0..1)
 233 // contributed by Matthias Troyer
 234 // for p=55, q=24 originally by G. J. Mitchell and D. P. Moore 1958
 235
 236 /**
 237  * Instantiations of class template @c lagged_fibonacci_01 model a
 238  * \pseudo_random_number_generator. It uses a lagged Fibonacci
 239  * algorithm with two lags @c p and @c q, evaluated in floating-point
 240  * arithmetic: x(i) = x(i-p) + x(i-q) (mod 1) with p > q. See
 241  *
 242  *  @blockquote
 243  *  "Uniform random number generators for supercomputers", Richard Brent,
 244  *  Proc. of Fifth Australian Supercomputer Conference, Melbourne,
 245  *  Dec. 1992, pp. 704-706.
 246  *  @endblockquote
 247  *
 248  * @xmlnote
 249  * The quality of the generator crucially depends on the choice
 250  * of the parameters. User code should employ one of the sensibly
 251  * parameterized generators such as \lagged_fibonacci607 instead.
 252  * @endxmlnote
 253  *
 254  * The generator requires considerable amounts of memory for the storage
 255  * of its state array. For example, \lagged_fibonacci607 requires about
 256  * 4856 bytes and \lagged_fibonacci44497 requires about 350 KBytes.
 257  */
 258 template<class RealType, int w, unsigned int p, unsigned int q>
 259 class lagged_fibonacci_01_engine
 260 {
 261 public:
 262     typedef RealType result_type;
 263     BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
 264     BOOST_STATIC_CONSTANT(int, word_size = w);
 265     BOOST_STATIC_CONSTANT(unsigned int, long_lag = p);
 266     BOOST_STATIC_CONSTANT(unsigned int, short_lag = q);
 267
 268     BOOST_STATIC_CONSTANT(boost::uint32_t, default_seed = 331u);
 269
 270     /** Constructs a @c lagged_fibonacci_01 generator and calls @c seed(). */
 271     lagged_fibonacci_01_engine() { seed(); }
 272     /** Constructs a @c lagged_fibonacci_01 generator and calls @c seed(value). */
 273     BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(lagged_fibonacci_01_engine, uint32_t, value)
 274     { seed(value); }
 275     /** Constructs a @c lagged_fibonacci_01 generator and calls @c seed(gen). */
 276     BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(lagged_fibonacci_01_engine, SeedSeq, seq)
 277     { seed(seq); }
 278     template<class It> lagged_fibonacci_01_engine(It& first, It last)
 279     { seed(first, last); }
 280
 281     // compiler-generated copy ctor and assignment operator are fine
 282
 283     /** Calls seed(default_seed). */
 284     void seed() { seed(default_seed); }
 285
 286     /**
 287      * Constructs a \minstd_rand0 generator with the constructor parameter
 288      * value and calls seed with it. Distinct seeds in the range
 289      * [1, 2147483647) will produce generators with different states. Other
 290      * seeds will be equivalent to some seed within this range. See
 291      * \linear_congruential_engine for details.
 292      */
 293     BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(lagged_fibonacci_01_engine, boost::uint32_t, value)
 294     {
 295         minstd_rand0 intgen(value);
 296         detail::generator_seed_seq<minstd_rand0> gen(intgen);
 297         seed(gen);
 298     }
 299
 300     /**
 301      * Seeds this @c lagged_fibonacci_01_engine using values produced by
 302      * @c seq.generate.
 303      */
 304     BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(lagged_fibonacci_01_engine, SeedSeq, seq)
 305     {
 306         detail::seed_array_real<w>(seq, x);
 307         i = long_lag;
 308     }
 309
 310     /**
 311      * Seeds this @c lagged_fibonacci_01_engine using values from the
 312      * iterator range [first, last).  If there are not enough elements
 313      * in the range, throws @c std::invalid_argument.
 314      */
 315     template<class It>
 316     void seed(It& first, It last)
 317     {
 318         detail::fill_array_real<w>(first, last, x);
 319         i = long_lag;
 320     }
 321
 322     /** Returns the smallest value that the generator can produce. */
 323     static result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () { return result_type(0); }
 324     /** Returns the upper bound of the generators outputs. */
 325     static result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () { return result_type(1); }
 326
 327     /** Returns the next value of the generator. */
 328     result_type operator()()
 329     {
 330         if(i >= long_lag)
 331             fill();
 332         return x[i++];
 333     }
 334
 335     /** Fills a range with random values */
 336     template<class Iter>
 337     void generate(Iter first, Iter last)
 338     { return detail::generate_from_real(*this, first, last); }
 339
 340     /** Advances the state of the generator by @c z. */
 341     void discard(boost::uintmax_t z)
 342     {
 343         for(boost::uintmax_t j = 0; j < z; ++j) {
 344             (*this)();
 345         }
 346     }
 347
 348     /**
 349      * Writes the textual representation of the generator to a @c std::ostream.
 350      */
 351     BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, lagged_fibonacci_01_engine, f)
 352     {
 353         // allow for Koenig lookup
 354         using std::pow;
 355         os << f.i;
 356         std::ios_base::fmtflags oldflags = os.flags(os.dec | os.fixed | os.left);
 357         for(unsigned int i = 0; i < f.long_lag; ++i)
 358             os << ' ' << f.x[i] * f.modulus();
 359         os.flags(oldflags);
 360         return os;
 361     }
 362
 363     /**
 364      * Reads the textual representation of the generator from a @c std::istream.
 365      */
 366     BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, lagged_fibonacci_01_engine, f)
 367     {
 368         is >> f.i;
 369         for(unsigned int i = 0; i < f.long_lag; ++i) {
 370             typename lagged_fibonacci_01_engine::result_type value;
 371             is >> std::ws >> value;
 372             f.x[i] = value / f.modulus();
 373         }
 374         return is;
 375     }
 376
 377     /**
 378      * Returns true if the two generators will produce identical
 379      * sequences of outputs.
 380      */
 381     BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(lagged_fibonacci_01_engine, x, y)
 382     { return x.i == y.i && std::equal(x.x, x.x+long_lag, y.x); }
 383
 384     /**
 385      * Returns true if the two generators will produce different
 386      * sequences of outputs.
 387      */
 388     BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(lagged_fibonacci_01_engine)
 389
 390 private:
 391     /// \cond show_private
 392     void fill();
 393     static RealType modulus()
 394     {
 395         using std::pow;
 396         return pow(RealType(2), word_size);
 397     }
 398     /// \endcond
 399     unsigned int i;
 400     RealType x[long_lag];
 401 };
 402
 403 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
 404 //  A definition is required even for integral static constants
 405 template<class RealType, int w, unsigned int p, unsigned int q>
 406 const bool lagged_fibonacci_01_engine<RealType, w, p, q>::has_fixed_range;
 407 template<class RealType, int w, unsigned int p, unsigned int q>
 408 const unsigned int lagged_fibonacci_01_engine<RealType, w, p, q>::long_lag;
 409 template<class RealType, int w, unsigned int p, unsigned int q>
 410 const unsigned int lagged_fibonacci_01_engine<RealType, w, p, q>::short_lag;
 411 template<class RealType, int w, unsigned int p, unsigned int q>
 412 const int lagged_fibonacci_01_engine<RealType,w,p,q>::word_size;
 413 template<class RealType, int w, unsigned int p, unsigned int q>
 414 const boost::uint32_t lagged_fibonacci_01_engine<RealType,w,p,q>::default_seed;
 415 #endif
 416
 417 /// \cond show_private
 418 template<class RealType, int w, unsigned int p, unsigned int q>
 419 void lagged_fibonacci_01_engine<RealType, w, p, q>::fill()
 420 {
 421     // two loops to avoid costly modulo operations
 422     {  // extra scope for MSVC brokenness w.r.t. for scope
 423     for(unsigned int j = 0; j < short_lag; ++j) {
 424         RealType t = x[j] + x[j+(long_lag-short_lag)];
 425         if(t >= RealType(1))
 426             t -= RealType(1);
 427         x[j] = t;
 428     }
 429     }
 430     for(unsigned int j = short_lag; j < long_lag; ++j) {
 431         RealType t = x[j] + x[j-short_lag];
 432         if(t >= RealType(1))
 433             t -= RealType(1);
 434         x[j] = t;
 435     }
 436     i = 0;
 437 }
 438 /// \endcond
 439
 440 /// \cond show_deprecated
 441
 442 // provided for backwards compatibility
 443 template<class RealType, int w, unsigned int p, unsigned int q>
 444 class lagged_fibonacci_01 : public lagged_fibonacci_01_engine<RealType, w, p, q>
 445 {
 446     typedef lagged_fibonacci_01_engine<RealType, w, p, q> base_type;
 447 public:
 448     lagged_fibonacci_01() {}
 449     BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(lagged_fibonacci_01, boost::uint32_t, val)
 450     { this->seed(val); }
 451     BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(lagged_fibonacci_01, SeedSeq, seq)
 452     { this->seed(seq); }
 453     template<class It>
 454     lagged_fibonacci_01(It& first, It last) : base_type(first, last) {}
 455 };
 456
 457 /// \endcond
 458
 459 namespace detail {
 460
 461 template<class Engine>
 462 struct generator_bits;
 463
 464 template<class RealType, int w, unsigned int p, unsigned int q>
 465 struct generator_bits<lagged_fibonacci_01_engine<RealType, w, p, q> >
 466 {
 467     static std::size_t value() { return w; }
 468 };
 469
 470 template<class RealType, int w, unsigned int p, unsigned int q>
 471 struct generator_bits<lagged_fibonacci_01<RealType, w, p, q> >
 472 {
 473     static std::size_t value() { return w; }
 474 };
 475
 476 }
 477
 478 #ifdef BOOST_RANDOM_DOXYGEN
 479 namespace detail {
 480 /**
 481  * The specializations lagged_fibonacci607 ... lagged_fibonacci44497
 482  * use well tested lags.
 483  *
 484  * See
 485  *
 486  *  @blockquote
 487  *  "On the Periods of Generalized Fibonacci Recurrences", Richard P. Brent
 488  *  Computer Sciences Laboratory Australian National University, December 1992
 489  *  @endblockquote
 490  *
 491  * The lags used here can be found in
 492  *
 493  *  @blockquote
 494  *  "Uniform random number generators for supercomputers", Richard Brent,
 495  *  Proc. of Fifth Australian Supercomputer Conference, Melbourne,
 496  *  Dec. 1992, pp. 704-706.
 497  *  @endblockquote
 498  */
 499 struct lagged_fibonacci_doc {};
 500 }
 501 #endif
 502
 503 /** @copydoc boost::random::detail::lagged_fibonacci_doc */
 504 typedef lagged_fibonacci_01_engine<double, 48, 607, 273> lagged_fibonacci607;
 505 /** @copydoc boost::random::detail::lagged_fibonacci_doc */
 506 typedef lagged_fibonacci_01_engine<double, 48, 1279, 418> lagged_fibonacci1279;
 507 /** @copydoc boost::random::detail::lagged_fibonacci_doc */
 508 typedef lagged_fibonacci_01_engine<double, 48, 2281, 1252> lagged_fibonacci2281;
 509 /** @copydoc boost::random::detail::lagged_fibonacci_doc */
 510 typedef lagged_fibonacci_01_engine<double, 48, 3217, 576> lagged_fibonacci3217;
 511 /** @copydoc boost::random::detail::lagged_fibonacci_doc */
 512 typedef lagged_fibonacci_01_engine<double, 48, 4423, 2098> lagged_fibonacci4423;
 513 /** @copydoc boost::random::detail::lagged_fibonacci_doc */
 514 typedef lagged_fibonacci_01_engine<double, 48, 9689, 5502> lagged_fibonacci9689;
 515 /** @copydoc boost::random::detail::lagged_fibonacci_doc */
 516 typedef lagged_fibonacci_01_engine<double, 48, 19937, 9842> lagged_fibonacci19937;
 517 /** @copydoc boost::random::detail::lagged_fibonacci_doc */
 518 typedef lagged_fibonacci_01_engine<double, 48, 23209, 13470> lagged_fibonacci23209;
 519 /** @copydoc boost::random::detail::lagged_fibonacci_doc */
 520 typedef lagged_fibonacci_01_engine<double, 48, 44497, 21034> lagged_fibonacci44497;
 521
 522 } // namespace random
 523
 524 using random::lagged_fibonacci607;
 525 using random::lagged_fibonacci1279;
 526 using random::lagged_fibonacci2281;
 527 using random::lagged_fibonacci3217;
 528 using random::lagged_fibonacci4423;
 529 using random::lagged_fibonacci9689;
 530 using random::lagged_fibonacci19937;
 531 using random::lagged_fibonacci23209;
 532 using random::lagged_fibonacci44497;
 533
 534 } // namespace boost
 535
 536 #endif // BOOST_RANDOM_LAGGED_FIBONACCI_HPP