1 /* boost random/lagged_fibonacci.hpp header file
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)
8 * See http://www.boost.org for most recent version including documentation.
10 * $Id: lagged_fibonacci.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $
13 * 2001-02-18 moved to individual header files
16 #ifndef BOOST_RANDOM_LAGGED_FIBONACCI_HPP
17 #define BOOST_RANDOM_LAGGED_FIBONACCI_HPP
19 #include <boost/config/no_tr1/cmath.hpp>
21 #include <algorithm> // std::max
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/detail/workaround.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/pass_through_engine.hpp>
37 #if BOOST_WORKAROUND(_MSC_FULL_VER, BOOST_TESTED_AT(13102292)) && BOOST_MSVC > 1300
38 # define BOOST_RANDOM_EXTRACT_LF
41 #if defined(__APPLE_CC__) && defined(__GNUC__) && (__GNUC__ == 3) && (__GNUC_MINOR__ <= 3)
42 # define BOOST_RANDOM_EXTRACT_LF
45 # ifdef BOOST_RANDOM_EXTRACT_LF
48 template<class IStream, class F, class RealType>
50 extract_lagged_fibonacci_01(
58 for(unsigned int i = 0; i < f.long_lag; ++i)
61 is >> value >> std::ws;
62 x[i] = value / modulus;
67 template<class IStream, class F, class UIntType>
69 extract_lagged_fibonacci(
76 for(unsigned int i = 0; i < f.long_lag; ++i)
77 is >> x[i] >> std::ws;
84 * Instantiations of class template \lagged_fibonacci model a
85 * \pseudo_random_number_generator. It uses a lagged Fibonacci
86 * algorithm with two lags @c p and @c q:
87 * x(i) = x(i-p) + x(i-q) (mod 2<sup>w</sup>) with p > q.
89 template<class UIntType, int w, unsigned int p, unsigned int q,
91 class lagged_fibonacci
94 typedef UIntType result_type;
95 BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
96 BOOST_STATIC_CONSTANT(int, word_size = w);
97 BOOST_STATIC_CONSTANT(unsigned int, long_lag = p);
98 BOOST_STATIC_CONSTANT(unsigned int, short_lag = q);
101 * Returns: the smallest value that the generator can produce
103 result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; }
105 * Returns: the largest value that the generator can produce
107 result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return wordmask; }
110 * Creates a new @c lagged_fibonacci generator and calls @c seed()
112 lagged_fibonacci() { init_wordmask(); seed(); }
114 * Creates a new @c lagged_fibonacci generator and calls @c seed(value)
116 explicit lagged_fibonacci(uint32_t value) { init_wordmask(); seed(value); }
118 * Creates a new @c lagged_fibonacci generator and calls @c seed(first, last)
120 template<class It> lagged_fibonacci(It& first, It last)
121 { init_wordmask(); seed(first, last); }
122 // compiler-generated copy ctor and assignment operator are fine
125 /// \cond hide_private_members
129 for(int j = 0; j < w; ++j)
130 wordmask |= (1u << j);
136 * Sets the state of the generator to values produced by
137 * a \minstd_rand generator.
139 void seed(uint32_t value = 331u)
141 minstd_rand0 gen(value);
142 for(unsigned int j = 0; j < long_lag; ++j)
143 x[j] = gen() & wordmask;
148 * Sets the state of the generator to values from the iterator
149 * range [first, last). If there are not enough elements in the
150 * range [first, last) throws @c std::invalid_argument.
153 void seed(It& first, It last)
155 // word size could be smaller than the seed values
157 for(j = 0; j < long_lag && first != last; ++j, ++first)
158 x[j] = *first & wordmask;
160 if(first == last && j < long_lag)
161 throw std::invalid_argument("lagged_fibonacci::seed");
165 * Returns: the next value of the generator
167 result_type operator()()
174 static bool validation(result_type x)
179 #ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
181 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
182 template<class CharT, class Traits>
183 friend std::basic_ostream<CharT,Traits>&
184 operator<<(std::basic_ostream<CharT,Traits>& os, const lagged_fibonacci& f)
187 for(unsigned int i = 0; i < f.long_lag; ++i)
192 template<class CharT, class Traits>
193 friend std::basic_istream<CharT, Traits>&
194 operator>>(std::basic_istream<CharT, Traits>& is, lagged_fibonacci& f)
196 # ifdef BOOST_RANDOM_EXTRACT_LF
197 return detail::extract_lagged_fibonacci(is, f, f.i, f.x);
199 is >> f.i >> std::ws;
200 for(unsigned int i = 0; i < f.long_lag; ++i)
201 is >> f.x[i] >> std::ws;
207 friend bool operator==(const lagged_fibonacci& x, const lagged_fibonacci& y)
208 { return x.i == y.i && std::equal(x.x, x.x+long_lag, y.x); }
209 friend bool operator!=(const lagged_fibonacci& x,
210 const lagged_fibonacci& y)
211 { return !(x == y); }
213 // Use a member function; Streamable concept not supported.
214 bool operator==(const lagged_fibonacci& rhs) const
215 { return i == rhs.i && std::equal(x, x+long_lag, rhs.x); }
216 bool operator!=(const lagged_fibonacci& rhs) const
217 { return !(*this == rhs); }
221 /// \cond hide_private_members
227 UIntType x[long_lag];
230 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
231 // A definition is required even for integral static constants
232 template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
233 const bool lagged_fibonacci<UIntType, w, p, q, val>::has_fixed_range;
234 template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
235 const unsigned int lagged_fibonacci<UIntType, w, p, q, val>::long_lag;
236 template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
237 const unsigned int lagged_fibonacci<UIntType, w, p, q, val>::short_lag;
240 /// \cond hide_private_members
242 template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
243 void lagged_fibonacci<UIntType, w, p, q, val>::fill()
245 // two loops to avoid costly modulo operations
246 { // extra scope for MSVC brokenness w.r.t. for scope
247 for(unsigned int j = 0; j < short_lag; ++j)
248 x[j] = (x[j] + x[j+(long_lag-short_lag)]) & wordmask;
250 for(unsigned int j = short_lag; j < long_lag; ++j)
251 x[j] = (x[j] + x[j-short_lag]) & wordmask;
257 // lagged Fibonacci generator for the range [0..1)
258 // contributed by Matthias Troyer
259 // for p=55, q=24 originally by G. J. Mitchell and D. P. Moore 1958
261 template<class T, unsigned int p, unsigned int q>
262 struct fibonacci_validation
264 BOOST_STATIC_CONSTANT(bool, is_specialized = false);
265 static T value() { return 0; }
266 static T tolerance() { return 0; }
269 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
270 // A definition is required even for integral static constants
271 template<class T, unsigned int p, unsigned int q>
272 const bool fibonacci_validation<T, p, q>::is_specialized;
275 #define BOOST_RANDOM_FIBONACCI_VAL(T,P,Q,V,E) \
277 struct fibonacci_validation<T, P, Q> \
279 BOOST_STATIC_CONSTANT(bool, is_specialized = true); \
280 static T value() { return V; } \
281 static T tolerance() \
282 { return (std::max)(E, static_cast<T>(5*std::numeric_limits<T>::epsilon())); } \
284 // (The extra static_cast<T> in the std::max call above is actually
285 // unnecessary except for HP aCC 1.30, which claims that
286 // numeric_limits<double>::epsilon() doesn't actually return a double.)
288 BOOST_RANDOM_FIBONACCI_VAL(double, 607, 273, 0.4293817707235914, 1e-14)
289 BOOST_RANDOM_FIBONACCI_VAL(double, 1279, 418, 0.9421630240437659, 1e-14)
290 BOOST_RANDOM_FIBONACCI_VAL(double, 2281, 1252, 0.1768114046909004, 1e-14)
291 BOOST_RANDOM_FIBONACCI_VAL(double, 3217, 576, 0.1956232694868209, 1e-14)
292 BOOST_RANDOM_FIBONACCI_VAL(double, 4423, 2098, 0.9499762202147172, 1e-14)
293 BOOST_RANDOM_FIBONACCI_VAL(double, 9689, 5502, 0.05737836943695162, 1e-14)
294 BOOST_RANDOM_FIBONACCI_VAL(double, 19937, 9842, 0.5076528587449834, 1e-14)
295 BOOST_RANDOM_FIBONACCI_VAL(double, 23209, 13470, 0.5414473810619185, 1e-14)
296 BOOST_RANDOM_FIBONACCI_VAL(double, 44497,21034, 0.254135073399297, 1e-14)
298 #undef BOOST_RANDOM_FIBONACCI_VAL
303 * Instantiations of class template @c lagged_fibonacci_01 model a
304 * \pseudo_random_number_generator. It uses a lagged Fibonacci
305 * algorithm with two lags @c p and @c q, evaluated in floating-point
306 * arithmetic: x(i) = x(i-p) + x(i-q) (mod 1) with p > q. See
309 * "Uniform random number generators for supercomputers", Richard Brent,
310 * Proc. of Fifth Australian Supercomputer Conference, Melbourne,
311 * Dec. 1992, pp. 704-706.
315 * The quality of the generator crucially depends on the choice
316 * of the parameters. User code should employ one of the sensibly
317 * parameterized generators such as \lagged_fibonacci607 instead.
320 * The generator requires considerable amounts of memory for the storage
321 * of its state array. For example, \lagged_fibonacci607 requires about
322 * 4856 bytes and \lagged_fibonacci44497 requires about 350 KBytes.
324 template<class RealType, int w, unsigned int p, unsigned int q>
325 class lagged_fibonacci_01
328 typedef RealType result_type;
329 BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
330 BOOST_STATIC_CONSTANT(int, word_size = w);
331 BOOST_STATIC_CONSTANT(unsigned int, long_lag = p);
332 BOOST_STATIC_CONSTANT(unsigned int, short_lag = q);
334 /** Constructs a @c lagged_fibonacci_01 generator and calls @c seed(). */
335 lagged_fibonacci_01() { init_modulus(); seed(); }
336 /** Constructs a @c lagged_fibonacci_01 generator and calls @c seed(value). */
337 BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(lagged_fibonacci_01, uint32_t, value)
338 { init_modulus(); seed(value); }
339 /** Constructs a @c lagged_fibonacci_01 generator and calls @c seed(gen). */
340 BOOST_RANDOM_DETAIL_GENERATOR_CONSTRUCTOR(lagged_fibonacci_01, Generator, gen)
341 { init_modulus(); seed(gen); }
342 template<class It> lagged_fibonacci_01(It& first, It last)
343 { init_modulus(); seed(first, last); }
344 // compiler-generated copy ctor and assignment operator are fine
347 /// \cond hide_private_members
350 #ifndef BOOST_NO_STDC_NAMESPACE
351 // allow for Koenig lookup
354 _modulus = pow(RealType(2), word_size);
359 /** Calls seed(331u). */
360 void seed() { seed(331u); }
362 * Constructs a \minstd_rand0 generator with the constructor parameter
363 * value and calls seed with it. Distinct seeds in the range
364 * [1, 2147483647) will produce generators with different states. Other
365 * seeds will be equivalent to some seed within this range. See
366 * \linear_congruential for details.
368 BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(lagged_fibonacci_01, uint32_t, value)
370 minstd_rand0 intgen(value);
375 * Sets the state of this @c lagged_fibonacci_01 to the values returned
376 * by p invocations of \uniform_01<code>\<RealType\>()(gen)</code>.
378 * Complexity: Exactly p invocations of gen.
380 BOOST_RANDOM_DETAIL_GENERATOR_SEED(lagged_fibonacci, Generator, gen)
382 // use pass-by-reference, but wrap argument in pass_through_engine
383 typedef detail::pass_through_engine<Generator&> ref_gen;
384 uniform_01<ref_gen, RealType> gen01 =
385 uniform_01<ref_gen, RealType>(ref_gen(gen));
386 // I could have used std::generate_n, but it takes "gen" by value
387 for(unsigned int j = 0; j < long_lag; ++j)
393 void seed(It& first, It last)
395 #ifndef BOOST_NO_STDC_NAMESPACE
396 // allow for Koenig lookup
400 unsigned long mask = ~((~0u) << (w%32)); // now lowest w bits set
401 RealType two32 = pow(RealType(2), 32);
403 for(j = 0; j < long_lag && first != last; ++j) {
405 for(int k = 0; k < w/32 && first != last; ++k, ++first)
406 x[j] += *first / pow(two32,k+1);
407 if(first != last && mask != 0) {
408 x[j] += fmod((*first & mask) / _modulus, RealType(1));
413 if(first == last && j < long_lag)
414 throw std::invalid_argument("lagged_fibonacci_01::seed");
417 result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(0); }
418 result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(1); }
420 result_type operator()()
427 static bool validation(result_type x)
429 result_type v = fibonacci_validation<result_type, p, q>::value();
430 result_type epsilon = fibonacci_validation<result_type, p, q>::tolerance();
431 // std::abs is a source of trouble: sometimes, it's not overloaded
432 // for double, plus the usual namespace std noncompliance -> avoid it
434 // return abs(x - v) < 5 * epsilon
435 return x > v - epsilon && x < v + epsilon;
438 #ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
440 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
441 template<class CharT, class Traits>
442 friend std::basic_ostream<CharT,Traits>&
443 operator<<(std::basic_ostream<CharT,Traits>& os, const lagged_fibonacci_01&f)
445 #ifndef BOOST_NO_STDC_NAMESPACE
446 // allow for Koenig lookup
450 std::ios_base::fmtflags oldflags = os.flags(os.dec | os.fixed | os.left);
451 for(unsigned int i = 0; i < f.long_lag; ++i)
452 os << f.x[i] * f._modulus << " ";
457 template<class CharT, class Traits>
458 friend std::basic_istream<CharT, Traits>&
459 operator>>(std::basic_istream<CharT, Traits>& is, lagged_fibonacci_01& f)
461 # ifdef BOOST_RANDOM_EXTRACT_LF
462 return detail::extract_lagged_fibonacci_01(is, f, f.i, f.x, f._modulus);
464 is >> f.i >> std::ws;
465 for(unsigned int i = 0; i < f.long_lag; ++i) {
466 typename lagged_fibonacci_01::result_type value;
467 is >> value >> std::ws;
468 f.x[i] = value / f._modulus;
475 friend bool operator==(const lagged_fibonacci_01& x,
476 const lagged_fibonacci_01& y)
477 { return x.i == y.i && std::equal(x.x, x.x+long_lag, y.x); }
478 friend bool operator!=(const lagged_fibonacci_01& x,
479 const lagged_fibonacci_01& y)
480 { return !(x == y); }
482 // Use a member function; Streamable concept not supported.
483 bool operator==(const lagged_fibonacci_01& rhs) const
484 { return i == rhs.i && std::equal(x, x+long_lag, rhs.x); }
485 bool operator!=(const lagged_fibonacci_01& rhs) const
486 { return !(*this == rhs); }
490 /// \cond hide_private_members
494 RealType x[long_lag];
498 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
499 // A definition is required even for integral static constants
500 template<class RealType, int w, unsigned int p, unsigned int q>
501 const bool lagged_fibonacci_01<RealType, w, p, q>::has_fixed_range;
502 template<class RealType, int w, unsigned int p, unsigned int q>
503 const unsigned int lagged_fibonacci_01<RealType, w, p, q>::long_lag;
504 template<class RealType, int w, unsigned int p, unsigned int q>
505 const unsigned int lagged_fibonacci_01<RealType, w, p, q>::short_lag;
506 template<class RealType, int w, unsigned int p, unsigned int q>
507 const int lagged_fibonacci_01<RealType,w,p,q>::word_size;
511 /// \cond hide_private_members
512 template<class RealType, int w, unsigned int p, unsigned int q>
513 void lagged_fibonacci_01<RealType, w, p, q>::fill()
515 // two loops to avoid costly modulo operations
516 { // extra scope for MSVC brokenness w.r.t. for scope
517 for(unsigned int j = 0; j < short_lag; ++j) {
518 RealType t = x[j] + x[j+(long_lag-short_lag)];
524 for(unsigned int j = short_lag; j < long_lag; ++j) {
525 RealType t = x[j] + x[j-short_lag];
534 } // namespace random
536 #ifdef BOOST_RANDOM_DOXYGEN
539 * The specializations lagged_fibonacci607 ... lagged_fibonacci44497
540 * use well tested lags.
545 * "On the Periods of Generalized Fibonacci Recurrences", Richard P. Brent
546 * Computer Sciences Laboratory Australian National University, December 1992
549 * The lags used here can be found in
552 * "Uniform random number generators for supercomputers", Richard Brent,
553 * Proc. of Fifth Australian Supercomputer Conference, Melbourne,
554 * Dec. 1992, pp. 704-706.
557 struct lagged_fibonacci_doc {};
562 * @copydoc boost::detail::lagged_fibonacci_doc
564 typedef random::lagged_fibonacci_01<double, 48, 607, 273> lagged_fibonacci607;
566 * @copydoc boost::detail::lagged_fibonacci_doc
568 typedef random::lagged_fibonacci_01<double, 48, 1279, 418> lagged_fibonacci1279;
570 * @copydoc boost::detail::lagged_fibonacci_doc
572 typedef random::lagged_fibonacci_01<double, 48, 2281, 1252> lagged_fibonacci2281;
574 * @copydoc boost::detail::lagged_fibonacci_doc
576 typedef random::lagged_fibonacci_01<double, 48, 3217, 576> lagged_fibonacci3217;
578 * @copydoc boost::detail::lagged_fibonacci_doc
580 typedef random::lagged_fibonacci_01<double, 48, 4423, 2098> lagged_fibonacci4423;
582 * @copydoc boost::detail::lagged_fibonacci_doc
584 typedef random::lagged_fibonacci_01<double, 48, 9689, 5502> lagged_fibonacci9689;
586 * @copydoc boost::detail::lagged_fibonacci_doc
588 typedef random::lagged_fibonacci_01<double, 48, 19937, 9842> lagged_fibonacci19937;
590 * @copydoc boost::detail::lagged_fibonacci_doc
592 typedef random::lagged_fibonacci_01<double, 48, 23209, 13470> lagged_fibonacci23209;
594 * @copydoc boost::detail::lagged_fibonacci_doc
596 typedef random::lagged_fibonacci_01<double, 48, 44497, 21034> lagged_fibonacci44497;
599 // It is possible to partially specialize uniform_01<> on lagged_fibonacci_01<>
600 // to help the compiler generate efficient code. For GCC, this seems useless,
601 // because GCC optimizes (x-0)/(1-0) to (x-0). This is good enough for now.
605 #endif // BOOST_RANDOM_LAGGED_FIBONACCI_HPP