source: NonGTP/Boost/boost/random/lagged_fibonacci.hpp @ 857

Revision 857, 15.3 KB checked in by igarcia, 19 years ago (diff)
Line 
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,v 1.28 2005/05/21 15:57:00 dgregor Exp $
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 <cmath>
20#include <iostream>
21#include <algorithm>     // std::max
22#include <iterator>
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>
31
32namespace boost {
33namespace random {
34
35#if BOOST_WORKAROUND(_MSC_FULL_VER, BOOST_TESTED_AT(13102292)) && BOOST_MSVC > 1300
36#  define BOOST_RANDOM_EXTRACT_LF
37#endif
38
39#if defined(__APPLE_CC__) && defined(__GNUC__) && (__GNUC__ == 3) && (__GNUC_MINOR__ <= 3)
40#  define BOOST_RANDOM_EXTRACT_LF
41#endif
42
43#  ifdef BOOST_RANDOM_EXTRACT_LF
44namespace detail
45{
46  template<class IStream, class F, class RealType>
47  IStream&
48  extract_lagged_fibonacci_01(
49      IStream& is
50      , F const& f
51      , unsigned int& i
52      , RealType* x
53      , RealType modulus)
54  {
55        is >> i >> std::ws;
56        for(unsigned int i = 0; i < f.long_lag; ++i)
57        {
58            RealType value;
59            is >> value >> std::ws;
60            x[i] = value / modulus;
61        }
62        return is;
63  }
64
65  template<class IStream, class F, class UIntType>
66  IStream&
67  extract_lagged_fibonacci(
68      IStream& is
69      , F const& f
70      , unsigned int& i
71      , UIntType* x)
72  {
73      is >> i >> std::ws;
74      for(unsigned int i = 0; i < f.long_lag; ++i)
75          is >> x[i] >> std::ws;
76      return is;
77  }
78}
79#  endif
80
81template<class UIntType, int w, unsigned int p, unsigned int q,
82         UIntType val = 0>
83class lagged_fibonacci
84{
85public:
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);
91
92  result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; }
93  result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return wordmask; }
94
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
100
101private:
102  void init_wordmask()
103  {
104    wordmask = 0;
105    for(int i = 0; i < w; ++i)
106      wordmask |= (1u << i);
107  }
108
109public:
110  void seed(uint32_t value = 331u)
111  {
112    minstd_rand0 gen(value);
113    for(unsigned int j = 0; j < long_lag; ++j)
114      x[j] = gen() & wordmask;
115    i = long_lag;
116  }
117
118  template<class It>
119  void seed(It& first, It last)
120  {
121    // word size could be smaller than the seed values
122    unsigned int j;
123    for(j = 0; j < long_lag && first != last; ++j, ++first)
124      x[j] = *first & wordmask;
125    i = long_lag;
126    if(first == last && j < long_lag)
127      throw std::invalid_argument("lagged_fibonacci::seed");
128  }
129
130  result_type operator()()
131  {
132    if(i >= long_lag)
133      fill();
134    return x[i++];
135  }
136
137  static bool validation(result_type x)
138  {
139    return x == val;
140  }
141 
142#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
143
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)
148  {
149    os << f.i << " ";
150    for(unsigned int i = 0; i < f.long_lag; ++i)
151      os << f.x[i] << " ";
152    return os;
153  }
154
155  template<class CharT, class Traits>
156  friend std::basic_istream<CharT, Traits>&
157  operator>>(std::basic_istream<CharT, Traits>& is, lagged_fibonacci& f)
158  {
159# ifdef BOOST_RANDOM_EXTRACT_LF
160      return detail::extract_lagged_fibonacci(is, f, f.i, f.x);
161# else
162      is >> f.i >> std::ws;
163      for(unsigned int i = 0; i < f.long_lag; ++i)
164          is >> f.x[i] >> std::ws;
165      return is;
166# endif
167  }
168#endif
169
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); }
175#else
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); }
181#endif
182
183private:
184  void fill();
185  UIntType wordmask;
186  unsigned int i;
187  UIntType x[long_lag];
188};
189
190#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
191//  A definition is required even for integral static constants
192template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
193const bool lagged_fibonacci<UIntType, w, p, q, val>::has_fixed_range;
194template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
195const unsigned int lagged_fibonacci<UIntType, w, p, q, val>::long_lag;
196template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
197const unsigned int lagged_fibonacci<UIntType, w, p, q, val>::short_lag;
198#endif
199
200template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
201void lagged_fibonacci<UIntType, w, p, q, val>::fill()
202{
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;
207  }
208  for(unsigned int j = short_lag; j < long_lag; ++j)
209    x[j] = (x[j] + x[j-short_lag]) & wordmask;
210  i = 0;
211}
212
213
214
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
218
219template<class T, unsigned int p, unsigned int q>
220struct fibonacci_validation
221{
222  BOOST_STATIC_CONSTANT(bool, is_specialized = false);
223  static T value() { return 0; }
224  static T tolerance() { return 0; }
225};
226
227#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
228//  A definition is required even for integral static constants
229template<class T, unsigned int p, unsigned int q>
230const bool fibonacci_validation<T, p, q>::is_specialized;
231#endif
232
233#define BOOST_RANDOM_FIBONACCI_VAL(T,P,Q,V,E) \
234template<> \
235struct fibonacci_validation<T, P, Q>  \
236{                                     \
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())); } \
241};
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.)
245
246BOOST_RANDOM_FIBONACCI_VAL(double, 607, 273, 0.4293817707235914, 1e-14)
247BOOST_RANDOM_FIBONACCI_VAL(double, 1279, 418, 0.9421630240437659, 1e-14)
248BOOST_RANDOM_FIBONACCI_VAL(double, 2281, 1252, 0.1768114046909004, 1e-14)
249BOOST_RANDOM_FIBONACCI_VAL(double, 3217, 576, 0.1956232694868209, 1e-14)
250BOOST_RANDOM_FIBONACCI_VAL(double, 4423, 2098, 0.9499762202147172, 1e-14)
251BOOST_RANDOM_FIBONACCI_VAL(double, 9689, 5502, 0.05737836943695162, 1e-14)
252BOOST_RANDOM_FIBONACCI_VAL(double, 19937, 9842, 0.5076528587449834, 1e-14)
253BOOST_RANDOM_FIBONACCI_VAL(double, 23209, 13470, 0.5414473810619185, 1e-14)
254BOOST_RANDOM_FIBONACCI_VAL(double, 44497,21034, 0.254135073399297, 1e-14)
255
256#undef BOOST_RANDOM_FIBONACCI_VAL
257
258template<class RealType, int w, unsigned int p, unsigned int q>
259class lagged_fibonacci_01
260{
261public:
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  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
275
276private:
277  void init_modulus()
278  {
279#ifndef BOOST_NO_STDC_NAMESPACE
280    // allow for Koenig lookup
281    using std::pow;
282#endif
283    _modulus = pow(RealType(2), word_size);
284  }
285
286public:
287  void seed(uint32_t value = 331u)
288  {
289    minstd_rand0 intgen(value);
290    seed(intgen);
291  }
292
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)
298  {
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)
305      x[j] = gen01();
306    i = long_lag;
307  }
308
309  template<class It>
310  void seed(It& first, It last)
311  {
312#ifndef BOOST_NO_STDC_NAMESPACE
313    // allow for Koenig lookup
314    using std::fmod;
315    using std::pow;
316#endif
317    unsigned long mask = ~((~0u) << (w%32));   // now lowest w bits set
318    RealType two32 = pow(RealType(2), 32);
319    unsigned int j;
320    for(j = 0; j < long_lag && first != last; ++j, ++first) {
321      x[j] = RealType(0);
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));
326    }
327    i = long_lag;
328    if(first == last && j < long_lag)
329      throw std::invalid_argument("lagged_fibonacci_01::seed");
330  }
331
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); }
334
335  result_type operator()()
336  {
337    if(i >= long_lag)
338      fill();
339    return x[i++];
340  }
341
342  static bool validation(result_type x)
343  {
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
348    // using std::abs;
349    // return abs(x - v) < 5 * epsilon
350    return x > v - epsilon && x < v + epsilon;
351  }
352 
353#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
354
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)
359  {
360#ifndef BOOST_NO_STDC_NAMESPACE
361    // allow for Koenig lookup
362    using std::pow;
363#endif
364    os << f.i << " ";
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 << " ";
368    os.flags(oldflags);
369    return os;
370  }
371
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)
375    {
376# ifdef BOOST_RANDOM_EXTRACT_LF
377        return detail::extract_lagged_fibonacci_01(is, f, f.i, f.x, f._modulus);
378# else
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;
384        }
385        return is;
386# endif
387    }
388#endif
389
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); }
396#else
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); }
402#endif
403
404private:
405  void fill();
406  unsigned int i;
407  RealType x[long_lag];
408  RealType _modulus;
409};
410
411#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
412//  A definition is required even for integral static constants
413template<class RealType, int w, unsigned int p, unsigned int q>
414const bool lagged_fibonacci_01<RealType, w, p, q>::has_fixed_range;
415template<class RealType, int w, unsigned int p, unsigned int q>
416const unsigned int lagged_fibonacci_01<RealType, w, p, q>::long_lag;
417template<class RealType, int w, unsigned int p, unsigned int q>
418const unsigned int lagged_fibonacci_01<RealType, w, p, q>::short_lag;
419template<class RealType, int w, unsigned int p, unsigned int q>
420const int lagged_fibonacci_01<RealType,w,p,q>::word_size;
421
422#endif
423
424template<class RealType, int w, unsigned int p, unsigned int q>
425void lagged_fibonacci_01<RealType, w, p, q>::fill()
426{
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)];
431    if(t >= RealType(1))
432      t -= RealType(1);
433    x[j] = t;
434  }
435  }
436  for(unsigned int j = short_lag; j < long_lag; ++j) {
437    RealType t = x[j] + x[j-short_lag];
438    if(t >= RealType(1))
439      t -= RealType(1);
440    x[j] = t;
441  }
442  i = 0;
443}
444
445} // namespace random
446
447typedef random::lagged_fibonacci_01<double, 48, 607, 273> lagged_fibonacci607;
448typedef random::lagged_fibonacci_01<double, 48, 1279, 418> lagged_fibonacci1279;
449typedef random::lagged_fibonacci_01<double, 48, 2281, 1252> lagged_fibonacci2281;
450typedef random::lagged_fibonacci_01<double, 48, 3217, 576> lagged_fibonacci3217;
451typedef random::lagged_fibonacci_01<double, 48, 4423, 2098> lagged_fibonacci4423;
452typedef random::lagged_fibonacci_01<double, 48, 9689, 5502> lagged_fibonacci9689;
453typedef random::lagged_fibonacci_01<double, 48, 19937, 9842> lagged_fibonacci19937;
454typedef random::lagged_fibonacci_01<double, 48, 23209, 13470> lagged_fibonacci23209;
455typedef random::lagged_fibonacci_01<double, 48, 44497, 21034> lagged_fibonacci44497;
456
457
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.
461
462} // namespace boost
463
464#endif // BOOST_RANDOM_LAGGED_FIBONACCI_HPP
Note: See TracBrowser for help on using the repository browser.