← Back to branch summary

~ubuntu-branches/ubuntu/warty/aqsis/warty

« back to all changes in this revision

Viewing changes to boost/boost/random/lagged_fibonacci.hpp

  • Committer: Bazaar Package Importer
  • Author(s): LaMont Jones
  • Date: 2004-08-24 07:25:04 UTC
  • Revision ID: james.westby@ubuntu.com-20040824072504-zf993vnevvisdsvb
Tags: upstream-0.9.1
Import upstream version 0.9.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
boost/boost/random/lagged_fibonacci.hpp
 
1
/* boost random/lagged_fibonacci.hpp header file
 
2
 *
 
3
 * Copyright Jens Maurer 2000-2001
 
4
 * Permission to use, copy, modify, sell, and distribute this software
 
5
 * is hereby granted without fee provided that the above copyright notice
 
6
 * appears in all copies and that both that copyright notice and this
 
7
 * permission notice appear in supporting documentation,
 
8
 *
 
9
 * Jens Maurer makes no representations about the suitability of this
 
10
 * software for any purpose. It is provided "as is" without express or
 
11
 * implied warranty.
 
12
 *
 
13
 * See http://www.boost.org for most recent version including documentation.
 
14
 *
 
15
 * $Id: lagged_fibonacci.hpp,v 1.1 2004/02/27 03:16:46 pseudonym Exp $
 
16
 *
 
17
 * Revision history
 
18
 *  2001-02-18  moved to individual header files
 
19
 */
 
20
 
 
21
#ifndef BOOST_RANDOM_LAGGED_FIBONACCI_HPP
 
22
#define BOOST_RANDOM_LAGGED_FIBONACCI_HPP
 
23
 
 
24
#include <iostream>
 
25
#include <algorithm>     // std::max
 
26
#include <iterator>
 
27
#include <boost/config.hpp>
 
28
#include <boost/limits.hpp>
 
29
#include <boost/cstdint.hpp>
 
30
#include <boost/random/linear_congruential.hpp>
 
31
#include <boost/random/uniform_01.hpp>
 
32
 
 
33
namespace boost {
 
34
namespace random {
 
35
 
 
36
#  if BOOST_WORKAROUND(_MSC_FULL_VER, BOOST_TESTED_AT(13102292) && BOOST_MSVC > 1300)
 
37
namespace detail
 
38
{
 
39
  template<class IStream, class F, class RealType>
 
40
  IStream&
 
41
  extract_lagged_fibonacci_01(
 
42
      IStream& is
 
43
      , F const& f
 
44
      , unsigned int& i
 
45
      , RealType* x
 
46
      , RealType modulus)
 
47
  {
 
48
        is >> i >> std::ws;
 
49
        for(unsigned int i = 0; i < f.long_lag; ++i)
 
50
        {
 
51
            RealType value;
 
52
            is >> value >> std::ws;
 
53
            x[i] = value / modulus;
 
54
        }
 
55
        return is;
 
56
  }
 
57
 
 
58
  template<class IStream, class F, class UIntType>
 
59
  IStream&
 
60
  extract_lagged_fibonacci(
 
61
      IStream& is
 
62
      , F const& f
 
63
      , unsigned int& i
 
64
      , UIntType* x)
 
65
  {
 
66
      is >> i >> std::ws;
 
67
      for(unsigned int i = 0; i < f.long_lag; ++i)
 
68
          is >> x[i] >> std::ws;
 
69
      return is;
 
70
  }
 
71
}
 
72
#  endif
 
73
 
 
74
template<class UIntType, int w, unsigned int p, unsigned int q,
 
75
         UIntType val = 0>
 
76
class lagged_fibonacci
 
77
{
 
78
public:
 
79
  typedef UIntType result_type;
 
80
  BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
 
81
  BOOST_STATIC_CONSTANT(int, word_size = w);
 
82
  BOOST_STATIC_CONSTANT(unsigned int, long_lag = p);
 
83
  BOOST_STATIC_CONSTANT(unsigned int, short_lag = q);
 
84
 
 
85
  result_type min() const { return 0; }
 
86
  result_type max() const { return wordmask; }
 
87
 
 
88
  lagged_fibonacci() { init_wordmask(); seed(); }
 
89
  explicit lagged_fibonacci(uint32_t value) { init_wordmask(); seed(value); }
 
90
  template<class It> lagged_fibonacci(It& first, It last)
 
91
  { init_wordmask(); seed(first, last); }
 
92
  // compiler-generated copy ctor and assignment operator are fine
 
93
 
 
94
private:
 
95
  void init_wordmask()
 
96
  {
 
97
    wordmask = 0;
 
98
    for(int i = 0; i < w; ++i)
 
99
      wordmask |= (1u << i);
 
100
  }
 
101
 
 
102
public:
 
103
  void seed(uint32_t value = 331u)
 
104
  {
 
105
    minstd_rand0 gen(value);
 
106
    for(unsigned int j = 0; j < long_lag; ++j)
 
107
      x[j] = gen() & wordmask;
 
108
    i = long_lag;
 
109
  }
 
110
 
 
111
  template<class It>
 
112
  void seed(It& first, It last)
 
113
  {
 
114
    // word size could be smaller than the seed values
 
115
    unsigned int j;
 
116
    for(j = 0; j < long_lag && first != last; ++j, ++first)
 
117
      x[j] = *first & wordmask;
 
118
    i = long_lag;
 
119
    if(first == last && j < long_lag)
 
120
      throw std::invalid_argument("lagged_fibonacci::seed");
 
121
  }
 
122
 
 
123
  result_type operator()()
 
124
  {
 
125
    if(i >= long_lag)
 
126
      fill();
 
127
    return x[i++];
 
128
  }
 
129
 
 
130
  static bool validation(result_type x)
 
131
  {
 
132
    return x == val;
 
133
  }
 
134
  
 
135
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
 
136
 
 
137
#ifndef BOOST_NO_MEMBER_TEMPLATE_FRIENDS
 
138
  template<class CharT, class Traits>
 
139
  friend std::basic_ostream<CharT,Traits>&
 
140
  operator<<(std::basic_ostream<CharT,Traits>& os, const lagged_fibonacci& f)
 
141
  {
 
142
    os << f.i << " ";
 
143
    for(unsigned int i = 0; i < f.long_lag; ++i)
 
144
      os << f.x[i] << " ";
 
145
    return os;
 
146
  }
 
147
 
 
148
  template<class CharT, class Traits>
 
149
  friend std::basic_istream<CharT, Traits>&
 
150
  operator>>(std::basic_istream<CharT, Traits>& is, lagged_fibonacci& f)
 
151
  {
 
152
# if BOOST_WORKAROUND(_MSC_FULL_VER, BOOST_TESTED_AT(13102292)) && BOOST_MSVC > 1300
 
153
      return detail::extract_lagged_fibonacci(is, f, f.i, f.x);
 
154
# else
 
155
      is >> f.i >> std::ws;
 
156
      for(unsigned int i = 0; i < f.long_lag; ++i)
 
157
          is >> f.x[i] >> std::ws;
 
158
      return is;
 
159
# endif 
 
160
  }
 
161
#endif
 
162
 
 
163
  friend bool operator==(const lagged_fibonacci& x, const lagged_fibonacci& y)
 
164
  { return x.i == y.i && std::equal(x.x, x.x+long_lag, y.x); }
 
165
  friend bool operator!=(const lagged_fibonacci& x,
 
166
                         const lagged_fibonacci& y)
 
167
  { return !(x == y); }
 
168
#else
 
169
  // Use a member function; Streamable concept not supported.
 
170
  bool operator==(const lagged_fibonacci& rhs) const
 
171
  { return i == rhs.i && std::equal(x, x+long_lag, rhs.x); }
 
172
  bool operator!=(const lagged_fibonacci& rhs) const
 
173
  { return !(*this == rhs); }
 
174
#endif
 
175
 
 
176
private:
 
177
  void fill();
 
178
  UIntType wordmask;
 
179
  unsigned int i;
 
180
  UIntType x[long_lag];
 
181
};
 
182
 
 
183
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
 
184
//  A definition is required even for integral static constants
 
185
template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
 
186
const bool lagged_fibonacci<UIntType, w, p, q, val>::has_fixed_range;
 
187
template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
 
188
const unsigned int lagged_fibonacci<UIntType, w, p, q, val>::long_lag;
 
189
template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
 
190
const unsigned int lagged_fibonacci<UIntType, w, p, q, val>::short_lag;
 
191
#endif
 
192
 
 
193
template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
 
194
void lagged_fibonacci<UIntType, w, p, q, val>::fill()
 
195
{
 
196
  // two loops to avoid costly modulo operations
 
197
  {  // extra scope for MSVC brokenness w.r.t. for scope
 
198
  for(unsigned int j = 0; j < short_lag; ++j)
 
199
    x[j] = (x[j] + x[j+(long_lag-short_lag)]) & wordmask;
 
200
  }
 
201
  for(unsigned int j = short_lag; j < long_lag; ++j)
 
202
    x[j] = (x[j] + x[j-short_lag]) & wordmask;
 
203
  i = 0;
 
204
}
 
205
 
 
206
 
 
207
 
 
208
// lagged Fibonacci generator for the range [0..1)
 
209
// contributed by Matthias Troyer
 
210
// for p=55, q=24 originally by G. J. Mitchell and D. P. Moore 1958
 
211
 
 
212
template<class T, unsigned int p, unsigned int q>
 
213
struct fibonacci_validation
 
214
{
 
215
  BOOST_STATIC_CONSTANT(bool, is_specialized = false);
 
216
  static T value() { return 0; }
 
217
  static T tolerance() { return 0; }
 
218
};
 
219
 
 
220
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
 
221
//  A definition is required even for integral static constants
 
222
template<class T, unsigned int p, unsigned int q>
 
223
const bool fibonacci_validation<T, p, q>::is_specialized;
 
224
#endif
 
225
 
 
226
#define BOOST_RANDOM_FIBONACCI_VAL(T,P,Q,V,E) \
 
227
template<> \
 
228
struct fibonacci_validation<T, P, Q>  \
 
229
{                                     \
 
230
  BOOST_STATIC_CONSTANT(bool, is_specialized = true);     \
 
231
  static T value() { return V; }      \
 
232
  static T tolerance()                \
 
233
    { return std::max(E, static_cast<T>(5*std::numeric_limits<T>::epsilon())); } \
 
234
};
 
235
// (The extra static_cast<T> in the std::max call above is actually
 
236
// unnecessary except for HP aCC 1.30, which claims that
 
237
// numeric_limits<double>::epsilon() doesn't actually return a double.)
 
238
 
 
239
BOOST_RANDOM_FIBONACCI_VAL(double, 607, 273, 0.4293817707235914, 1e-14)
 
240
BOOST_RANDOM_FIBONACCI_VAL(double, 1279, 418, 0.9421630240437659, 1e-14)
 
241
BOOST_RANDOM_FIBONACCI_VAL(double, 2281, 1252, 0.1768114046909004, 1e-14)
 
242
BOOST_RANDOM_FIBONACCI_VAL(double, 3217, 576, 0.1956232694868209, 1e-14)
 
243
BOOST_RANDOM_FIBONACCI_VAL(double, 4423, 2098, 0.9499762202147172, 1e-14)
 
244
BOOST_RANDOM_FIBONACCI_VAL(double, 9689, 5502, 0.05737836943695162, 1e-14)
 
245
BOOST_RANDOM_FIBONACCI_VAL(double, 19937, 9842, 0.5076528587449834, 1e-14)
 
246
BOOST_RANDOM_FIBONACCI_VAL(double, 23209, 13470, 0.5414473810619185, 1e-14)
 
247
BOOST_RANDOM_FIBONACCI_VAL(double, 44497,21034, 0.254135073399297, 1e-14)
 
248
 
 
249
#undef BOOST_RANDOM_FIBONACCI_VAL
 
250
 
 
251
template<class RealType, int w, unsigned int p, unsigned int q>
 
252
class lagged_fibonacci_01
 
253
{
 
254
public:
 
255
  typedef RealType result_type;
 
256
  BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
 
257
  BOOST_STATIC_CONSTANT(int, word_size = w);
 
258
  BOOST_STATIC_CONSTANT(unsigned int, long_lag = p);
 
259
  BOOST_STATIC_CONSTANT(unsigned int, short_lag = q);
 
260
 
 
261
  lagged_fibonacci_01() { init_modulus(); seed(); }
 
262
  explicit lagged_fibonacci_01(uint32_t value) { init_modulus(); seed(value); }
 
263
  template<class Generator>
 
264
  explicit lagged_fibonacci_01(Generator & gen) { init_modulus(); seed(gen); }
 
265
  template<class It> lagged_fibonacci_01(It& first, It last)
 
266
  { init_modulus(); seed(first, last); }
 
267
  // compiler-generated copy ctor and assignment operator are fine
 
268
 
 
269
private:
 
270
  void init_modulus()
 
271
  {
 
272
#ifndef BOOST_NO_STDC_NAMESPACE
 
273
    // allow for Koenig lookup
 
274
    using std::pow;
 
275
#endif
 
276
    _modulus = pow(RealType(2), word_size);
 
277
  }
 
278
 
 
279
public:
 
280
  void seed(uint32_t value = 331u)
 
281
  {
 
282
    minstd_rand0 intgen(value);
 
283
    seed(intgen);
 
284
  }
 
285
 
 
286
  // For GCC, moving this function out-of-line prevents inlining, which may
 
287
  // reduce overall object code size.  However, MSVC does not grok
 
288
  // out-of-line template member functions.
 
289
  template<class Generator>
 
290
  void seed(Generator & gen)
 
291
  {
 
292
    uniform_01<Generator, RealType> gen01(gen);
 
293
    // I could have used std::generate_n, but it takes "gen" by value
 
294
    for(unsigned int j = 0; j < long_lag; ++j)
 
295
      x[j] = gen01();
 
296
    i = long_lag;
 
297
  }
 
298
 
 
299
  template<class It>
 
300
  void seed(It& first, It last)
 
301
  {
 
302
#ifndef BOOST_NO_STDC_NAMESPACE
 
303
    // allow for Koenig lookup
 
304
    using std::fmod;
 
305
    using std::pow;
 
306
#endif
 
307
    unsigned long mask = ~((~0u) << (w%32));   // now lowest w bits set
 
308
    RealType two32 = pow(RealType(2), 32);
 
309
    unsigned int j;
 
310
    for(j = 0; j < long_lag && first != last; ++j, ++first) {
 
311
      x[j] = RealType(0);
 
312
      for(int i = 0; i < w/32 && first != last; ++i, ++first)
 
313
        x[j] += *first / pow(two32,i+1);
 
314
      if(first != last && mask != 0)
 
315
        x[j] += fmod((*first & mask) / _modulus, RealType(1));
 
316
    }
 
317
    i = long_lag;
 
318
    if(first == last && j < long_lag)
 
319
      throw std::invalid_argument("lagged_fibonacci_01::seed");
 
320
  }
 
321
 
 
322
  result_type min() const { return result_type(0); }
 
323
  result_type max() const { return result_type(1); }
 
324
 
 
325
  result_type operator()()
 
326
  {
 
327
    if(i >= long_lag)
 
328
      fill();
 
329
    return x[i++];
 
330
  }
 
331
 
 
332
  static bool validation(result_type x)
 
333
  {
 
334
    result_type v = fibonacci_validation<result_type, p, q>::value();
 
335
    result_type epsilon = fibonacci_validation<result_type, p, q>::tolerance();
 
336
    // std::abs is a source of trouble: sometimes, it's not overloaded
 
337
    // for double, plus the usual namespace std noncompliance -> avoid it
 
338
    // using std::abs;
 
339
    // return abs(x - v) < 5 * epsilon
 
340
    return x > v - epsilon && x < v + epsilon;
 
341
  }
 
342
  
 
343
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
 
344
 
 
345
#ifndef BOOST_NO_MEMBER_TEMPLATE_FRIENDS
 
346
  template<class CharT, class Traits>
 
347
  friend std::basic_ostream<CharT,Traits>&
 
348
  operator<<(std::basic_ostream<CharT,Traits>& os, const lagged_fibonacci_01&f)
 
349
  {
 
350
#ifndef BOOST_NO_STDC_NAMESPACE
 
351
    // allow for Koenig lookup
 
352
    using std::pow;
 
353
#endif
 
354
    os << f.i << " ";
 
355
    std::ios_base::fmtflags oldflags = os.flags(os.dec | os.fixed | os.left); 
 
356
    for(unsigned int i = 0; i < f.long_lag; ++i)
 
357
      os << f.x[i] * f._modulus << " ";
 
358
    os.flags(oldflags);
 
359
    return os;
 
360
  }
 
361
 
 
362
  template<class CharT, class Traits>
 
363
  friend std::basic_istream<CharT, Traits>&
 
364
  operator>>(std::basic_istream<CharT, Traits>& is, lagged_fibonacci_01& f)
 
365
    {
 
366
# if BOOST_WORKAROUND(_MSC_FULL_VER, BOOST_TESTED_AT(13102292)) && BOOST_MSVC > 1300
 
367
        return detail::extract_lagged_fibonacci_01(is, f, f.i, f.x, f._modulus);
 
368
# else
 
369
        is >> f.i >> std::ws;
 
370
        for(unsigned int i = 0; i < f.long_lag; ++i) {
 
371
            RealType value;
 
372
            is >> value >> std::ws;
 
373
            f.x[i] = value / f._modulus;
 
374
        }
 
375
        return is;
 
376
# endif 
 
377
    }
 
378
#endif
 
379
 
 
380
  friend bool operator==(const lagged_fibonacci_01& x,
 
381
                         const lagged_fibonacci_01& y)
 
382
  { return x.i == y.i && std::equal(x.x, x.x+long_lag, y.x); }
 
383
  friend bool operator!=(const lagged_fibonacci_01& x,
 
384
                         const lagged_fibonacci_01& y)
 
385
  { return !(x == y); }
 
386
#else
 
387
  // Use a member function; Streamable concept not supported.
 
388
  bool operator==(const lagged_fibonacci_01& rhs) const
 
389
  { return i == rhs.i && std::equal(x, x+long_lag, rhs.x); }
 
390
  bool operator!=(const lagged_fibonacci_01& rhs) const
 
391
  { return !(*this == rhs); }
 
392
#endif
 
393
 
 
394
private:
 
395
  void fill();
 
396
  unsigned int i;
 
397
  RealType x[long_lag];
 
398
  RealType _modulus;
 
399
};
 
400
 
 
401
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
 
402
//  A definition is required even for integral static constants
 
403
template<class RealType, int w, unsigned int p, unsigned int q>
 
404
const bool lagged_fibonacci_01<RealType, w, p, q>::has_fixed_range;
 
405
template<class RealType, int w, unsigned int p, unsigned int q>
 
406
const unsigned int lagged_fibonacci_01<RealType, w, p, q>::long_lag;
 
407
template<class RealType, int w, unsigned int p, unsigned int q>
 
408
const unsigned int lagged_fibonacci_01<RealType, w, p, q>::short_lag;
 
409
 
 
410
#endif
 
411
 
 
412
template<class RealType, int w, unsigned int p, unsigned int q>
 
413
void lagged_fibonacci_01<RealType, w, p, q>::fill()
 
414
{
 
415
  // two loops to avoid costly modulo operations
 
416
  {  // extra scope for MSVC brokenness w.r.t. for scope
 
417
  for(unsigned int j = 0; j < short_lag; ++j) {
 
418
    RealType t = x[j] + x[j+(long_lag-short_lag)];
 
419
    if(t >= RealType(1))
 
420
      t -= RealType(1);
 
421
    x[j] = t;
 
422
  }
 
423
  }
 
424
  for(unsigned int j = short_lag; j < long_lag; ++j) {
 
425
    RealType t = x[j] + x[j-short_lag];
 
426
    if(t >= RealType(1))
 
427
      t -= RealType(1);
 
428
    x[j] = t;
 
429
  }
 
430
  i = 0;
 
431
}
 
432
 
 
433
} // namespace random
 
434
 
 
435
typedef random::lagged_fibonacci_01<double, 48, 607, 273> lagged_fibonacci607;
 
436
typedef random::lagged_fibonacci_01<double, 48, 1279, 418> lagged_fibonacci1279;
 
437
typedef random::lagged_fibonacci_01<double, 48, 2281, 1252> lagged_fibonacci2281;
 
438
typedef random::lagged_fibonacci_01<double, 48, 3217, 576> lagged_fibonacci3217;
 
439
typedef random::lagged_fibonacci_01<double, 48, 4423, 2098> lagged_fibonacci4423;
 
440
typedef random::lagged_fibonacci_01<double, 48, 9689, 5502> lagged_fibonacci9689;
 
441
typedef random::lagged_fibonacci_01<double, 48, 19937, 9842> lagged_fibonacci19937;
 
442
typedef random::lagged_fibonacci_01<double, 48, 23209, 13470> lagged_fibonacci23209;
 
443
typedef random::lagged_fibonacci_01<double, 48, 44497, 21034> lagged_fibonacci44497;
 
444
 
 
445
 
 
446
// It is possible to partially specialize uniform_01<> on lagged_fibonacci_01<>
 
447
// to help the compiler generate efficient code.  For GCC, this seems useless,
 
448
// because GCC optimizes (x-0)/(1-0) to (x-0).  This is good enough for now.
 
449
 
 
450
} // namespace boost
 
451
 
 
452
#endif // BOOST_RANDOM_LAGGED_FIBONACCI_HPP