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,v 1.28 2005/05/21 15:57:00 dgregor Exp $
13 * 2001-02-18 moved to individual header files
16 #ifndef BOOST_RANDOM_LAGGED_FIBONACCI_HPP
17 #define BOOST_RANDOM_LAGGED_FIBONACCI_HPP
21 #include <algorithm> // std::max
23 #include <cmath> // 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/pass_through_engine.hpp>
35 #if BOOST_WORKAROUND(_MSC_FULL_VER, BOOST_TESTED_AT(13102292)) && BOOST_MSVC > 1300
36 # define BOOST_RANDOM_EXTRACT_LF
39 #if defined(__APPLE_CC__) && defined(__GNUC__) && (__GNUC__ == 3) && (__GNUC_MINOR__ <= 3)
40 # define BOOST_RANDOM_EXTRACT_LF
43 # ifdef BOOST_RANDOM_EXTRACT_LF
46 template<class IStream, class F, class RealType>
48 extract_lagged_fibonacci_01(
56 for(unsigned int i = 0; i < f.long_lag; ++i)
59 is >> value >> std::ws;
60 x[i] = value / modulus;
65 template<class IStream, class F, class UIntType>
67 extract_lagged_fibonacci(
74 for(unsigned int i = 0; i < f.long_lag; ++i)
75 is >> x[i] >> std::ws;
81 template<class UIntType, int w, unsigned int p, unsigned int q,
83 class lagged_fibonacci
86 typedef UIntType result_type;
87 BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
88 BOOST_STATIC_CONSTANT(int, word_size = w);
89 BOOST_STATIC_CONSTANT(unsigned int, long_lag = p);
90 BOOST_STATIC_CONSTANT(unsigned int, short_lag = q);
92 result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; }
93 result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return wordmask; }
95 lagged_fibonacci() { init_wordmask(); seed(); }
96 explicit lagged_fibonacci(uint32_t value) { init_wordmask(); seed(value); }
97 template<class It> lagged_fibonacci(It& first, It last)
98 { init_wordmask(); seed(first, last); }
99 // compiler-generated copy ctor and assignment operator are fine
105 for(int i = 0; i < w; ++i)
106 wordmask |= (1u << i);
110 void seed(uint32_t value = 331u)
112 minstd_rand0 gen(value);
113 for(unsigned int j = 0; j < long_lag; ++j)
114 x[j] = gen() & wordmask;
119 void seed(It& first, It last)
121 // word size could be smaller than the seed values
123 for(j = 0; j < long_lag && first != last; ++j, ++first)
124 x[j] = *first & wordmask;
126 if(first == last && j < long_lag)
127 throw std::invalid_argument("lagged_fibonacci::seed");
130 result_type operator()()
137 static bool validation(result_type x)
142 #ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
144 #ifndef BOOST_NO_MEMBER_TEMPLATE_FRIENDS
145 template<class CharT, class Traits>
146 friend std::basic_ostream<CharT,Traits>&
147 operator<<(std::basic_ostream<CharT,Traits>& os, const lagged_fibonacci& f)
150 for(unsigned int i = 0; i < f.long_lag; ++i)
155 template<class CharT, class Traits>
156 friend std::basic_istream<CharT, Traits>&
157 operator>>(std::basic_istream<CharT, Traits>& is, lagged_fibonacci& f)
159 # ifdef BOOST_RANDOM_EXTRACT_LF
160 return detail::extract_lagged_fibonacci(is, f, f.i, f.x);
162 is >> f.i >> std::ws;
163 for(unsigned int i = 0; i < f.long_lag; ++i)
164 is >> f.x[i] >> std::ws;
170 friend bool operator==(const lagged_fibonacci& x, const lagged_fibonacci& y)
171 { return x.i == y.i && std::equal(x.x, x.x+long_lag, y.x); }
172 friend bool operator!=(const lagged_fibonacci& x,
173 const lagged_fibonacci& y)
174 { return !(x == y); }
176 // Use a member function; Streamable concept not supported.
177 bool operator==(const lagged_fibonacci& rhs) const
178 { return i == rhs.i && std::equal(x, x+long_lag, rhs.x); }
179 bool operator!=(const lagged_fibonacci& rhs) const
180 { return !(*this == rhs); }
187 UIntType x[long_lag];
190 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
191 // A definition is required even for integral static constants
192 template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
193 const bool lagged_fibonacci<UIntType, w, p, q, val>::has_fixed_range;
194 template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
195 const unsigned int lagged_fibonacci<UIntType, w, p, q, val>::long_lag;
196 template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
197 const unsigned int lagged_fibonacci<UIntType, w, p, q, val>::short_lag;
200 template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
201 void lagged_fibonacci<UIntType, w, p, q, val>::fill()
203 // two loops to avoid costly modulo operations
204 { // extra scope for MSVC brokenness w.r.t. for scope
205 for(unsigned int j = 0; j < short_lag; ++j)
206 x[j] = (x[j] + x[j+(long_lag-short_lag)]) & wordmask;
208 for(unsigned int j = short_lag; j < long_lag; ++j)
209 x[j] = (x[j] + x[j-short_lag]) & wordmask;
215 // lagged Fibonacci generator for the range [0..1)
216 // contributed by Matthias Troyer
217 // for p=55, q=24 originally by G. J. Mitchell and D. P. Moore 1958
219 template<class T, unsigned int p, unsigned int q>
220 struct fibonacci_validation
222 BOOST_STATIC_CONSTANT(bool, is_specialized = false);
223 static T value() { return 0; }
224 static T tolerance() { return 0; }
227 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
228 // A definition is required even for integral static constants
229 template<class T, unsigned int p, unsigned int q>
230 const bool fibonacci_validation<T, p, q>::is_specialized;
233 #define BOOST_RANDOM_FIBONACCI_VAL(T,P,Q,V,E) \
235 struct fibonacci_validation<T, P, Q> \
237 BOOST_STATIC_CONSTANT(bool, is_specialized = true); \
238 static T value() { return V; } \
239 static T tolerance() \
240 { return (std::max)(E, static_cast<T>(5*std::numeric_limits<T>::epsilon())); } \
242 // (The extra static_cast<T> in the std::max call above is actually
243 // unnecessary except for HP aCC 1.30, which claims that
244 // numeric_limits<double>::epsilon() doesn't actually return a double.)
246 BOOST_RANDOM_FIBONACCI_VAL(double, 607, 273, 0.4293817707235914, 1e-14)
247 BOOST_RANDOM_FIBONACCI_VAL(double, 1279, 418, 0.9421630240437659, 1e-14)
248 BOOST_RANDOM_FIBONACCI_VAL(double, 2281, 1252, 0.1768114046909004, 1e-14)
249 BOOST_RANDOM_FIBONACCI_VAL(double, 3217, 576, 0.1956232694868209, 1e-14)
250 BOOST_RANDOM_FIBONACCI_VAL(double, 4423, 2098, 0.9499762202147172, 1e-14)
251 BOOST_RANDOM_FIBONACCI_VAL(double, 9689, 5502, 0.05737836943695162, 1e-14)
252 BOOST_RANDOM_FIBONACCI_VAL(double, 19937, 9842, 0.5076528587449834, 1e-14)
253 BOOST_RANDOM_FIBONACCI_VAL(double, 23209, 13470, 0.5414473810619185, 1e-14)
254 BOOST_RANDOM_FIBONACCI_VAL(double, 44497,21034, 0.254135073399297, 1e-14)
256 #undef BOOST_RANDOM_FIBONACCI_VAL
258 template<class RealType, int w, unsigned int p, unsigned int q>
259 class lagged_fibonacci_01
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);
268 lagged_fibonacci_01() { init_modulus(); seed(); }
269 explicit lagged_fibonacci_01(uint32_t value) { init_modulus(); seed(value); }
270 template<class Generator>
271 explicit lagged_fibonacci_01(Generator & gen) { init_modulus(); seed(gen); }
272 template<class It> lagged_fibonacci_01(It& first, It last)
273 { init_modulus(); seed(first, last); }
274 // compiler-generated copy ctor and assignment operator are fine
279 #ifndef BOOST_NO_STDC_NAMESPACE
280 // allow for Koenig lookup
283 _modulus = pow(RealType(2), word_size);
287 void seed(uint32_t value = 331u)
289 minstd_rand0 intgen(value);
293 // For GCC, moving this function out-of-line prevents inlining, which may
294 // reduce overall object code size. However, MSVC does not grok
295 // out-of-line template member functions.
296 template<class Generator>
297 void seed(Generator & gen)
299 // use pass-by-reference, but wrap argument in pass_through_engine
300 typedef detail::pass_through_engine<Generator&> ref_gen;
301 uniform_01<ref_gen, RealType> gen01 =
302 uniform_01<ref_gen, RealType>(ref_gen(gen));
303 // I could have used std::generate_n, but it takes "gen" by value
304 for(unsigned int j = 0; j < long_lag; ++j)
310 void seed(It& first, It last)
312 #ifndef BOOST_NO_STDC_NAMESPACE
313 // allow for Koenig lookup
317 unsigned long mask = ~((~0u) << (w%32)); // now lowest w bits set
318 RealType two32 = pow(RealType(2), 32);
320 for(j = 0; j < long_lag && first != last; ++j, ++first) {
322 for(int i = 0; i < w/32 && first != last; ++i, ++first)
323 x[j] += *first / pow(two32,i+1);
324 if(first != last && mask != 0)
325 x[j] += fmod((*first & mask) / _modulus, RealType(1));
328 if(first == last && j < long_lag)
329 throw std::invalid_argument("lagged_fibonacci_01::seed");
332 result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(0); }
333 result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(1); }
335 result_type operator()()
342 static bool validation(result_type x)
344 result_type v = fibonacci_validation<result_type, p, q>::value();
345 result_type epsilon = fibonacci_validation<result_type, p, q>::tolerance();
346 // std::abs is a source of trouble: sometimes, it's not overloaded
347 // for double, plus the usual namespace std noncompliance -> avoid it
349 // return abs(x - v) < 5 * epsilon
350 return x > v - epsilon && x < v + epsilon;
353 #ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
355 #ifndef BOOST_NO_MEMBER_TEMPLATE_FRIENDS
356 template<class CharT, class Traits>
357 friend std::basic_ostream<CharT,Traits>&
358 operator<<(std::basic_ostream<CharT,Traits>& os, const lagged_fibonacci_01&f)
360 #ifndef BOOST_NO_STDC_NAMESPACE
361 // allow for Koenig lookup
365 std::ios_base::fmtflags oldflags = os.flags(os.dec | os.fixed | os.left);
366 for(unsigned int i = 0; i < f.long_lag; ++i)
367 os << f.x[i] * f._modulus << " ";
372 template<class CharT, class Traits>
373 friend std::basic_istream<CharT, Traits>&
374 operator>>(std::basic_istream<CharT, Traits>& is, lagged_fibonacci_01& f)
376 # ifdef BOOST_RANDOM_EXTRACT_LF
377 return detail::extract_lagged_fibonacci_01(is, f, f.i, f.x, f._modulus);
379 is >> f.i >> std::ws;
380 for(unsigned int i = 0; i < f.long_lag; ++i) {
381 typename lagged_fibonacci_01::result_type value;
382 is >> value >> std::ws;
383 f.x[i] = value / f._modulus;
390 friend bool operator==(const lagged_fibonacci_01& x,
391 const lagged_fibonacci_01& y)
392 { return x.i == y.i && std::equal(x.x, x.x+long_lag, y.x); }
393 friend bool operator!=(const lagged_fibonacci_01& x,
394 const lagged_fibonacci_01& y)
395 { return !(x == y); }
397 // Use a member function; Streamable concept not supported.
398 bool operator==(const lagged_fibonacci_01& rhs) const
399 { return i == rhs.i && std::equal(x, x+long_lag, rhs.x); }
400 bool operator!=(const lagged_fibonacci_01& rhs) const
401 { return !(*this == rhs); }
407 RealType x[long_lag];
411 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
412 // A definition is required even for integral static constants
413 template<class RealType, int w, unsigned int p, unsigned int q>
414 const bool lagged_fibonacci_01<RealType, w, p, q>::has_fixed_range;
415 template<class RealType, int w, unsigned int p, unsigned int q>
416 const unsigned int lagged_fibonacci_01<RealType, w, p, q>::long_lag;
417 template<class RealType, int w, unsigned int p, unsigned int q>
418 const unsigned int lagged_fibonacci_01<RealType, w, p, q>::short_lag;
419 template<class RealType, int w, unsigned int p, unsigned int q>
420 const int lagged_fibonacci_01<RealType,w,p,q>::word_size;
424 template<class RealType, int w, unsigned int p, unsigned int q>
425 void lagged_fibonacci_01<RealType, w, p, q>::fill()
427 // two loops to avoid costly modulo operations
428 { // extra scope for MSVC brokenness w.r.t. for scope
429 for(unsigned int j = 0; j < short_lag; ++j) {
430 RealType t = x[j] + x[j+(long_lag-short_lag)];
436 for(unsigned int j = short_lag; j < long_lag; ++j) {
437 RealType t = x[j] + x[j-short_lag];
445 } // namespace random
447 typedef random::lagged_fibonacci_01<double, 48, 607, 273> lagged_fibonacci607;
448 typedef random::lagged_fibonacci_01<double, 48, 1279, 418> lagged_fibonacci1279;
449 typedef random::lagged_fibonacci_01<double, 48, 2281, 1252> lagged_fibonacci2281;
450 typedef random::lagged_fibonacci_01<double, 48, 3217, 576> lagged_fibonacci3217;
451 typedef random::lagged_fibonacci_01<double, 48, 4423, 2098> lagged_fibonacci4423;
452 typedef random::lagged_fibonacci_01<double, 48, 9689, 5502> lagged_fibonacci9689;
453 typedef random::lagged_fibonacci_01<double, 48, 19937, 9842> lagged_fibonacci19937;
454 typedef random::lagged_fibonacci_01<double, 48, 23209, 13470> lagged_fibonacci23209;
455 typedef random::lagged_fibonacci_01<double, 48, 44497, 21034> lagged_fibonacci44497;
458 // It is possible to partially specialize uniform_01<> on lagged_fibonacci_01<>
459 // to help the compiler generate efficient code. For GCC, this seems useless,
460 // because GCC optimizes (x-0)/(1-0) to (x-0). This is good enough for now.
464 #endif // BOOST_RANDOM_LAGGED_FIBONACCI_HPP