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

Revision 857, 14.1 KB checked in by igarcia, 19 years ago (diff)
RevLine 
[857]1/* boost random/subtract_with_carry.hpp header file
2 *
3 * Copyright Jens Maurer 2002
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: subtract_with_carry.hpp,v 1.24 2005/05/21 15:57:00 dgregor Exp $
11 *
12 * Revision history
13 *  2002-03-02  created
14 */
15
16#ifndef BOOST_RANDOM_SUBTRACT_WITH_CARRY_HPP
17#define BOOST_RANDOM_SUBTRACT_WITH_CARRY_HPP
18
19#include <cmath>
20#include <iostream>
21#include <algorithm>     // std::equal
22#include <stdexcept>
23#include <cmath>         // std::pow
24#include <boost/config.hpp>
25#include <boost/limits.hpp>
26#include <boost/cstdint.hpp>
27#include <boost/static_assert.hpp>
28#include <boost/detail/workaround.hpp>
29#include <boost/random/linear_congruential.hpp>
30
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_SWC_01
37#endif
38
39#if defined(__APPLE_CC__) && defined(__GNUC__) && (__GNUC__ == 3) && (__GNUC_MINOR__ <= 3)
40#  define BOOST_RANDOM_EXTRACT_SWC_01
41#endif
42
43# ifdef BOOST_RANDOM_EXTRACT_SWC_01
44namespace detail
45{
46  template <class IStream, class SubtractWithCarry, class RealType>
47  void extract_subtract_with_carry_01(
48      IStream& is
49      , SubtractWithCarry& f
50      , RealType& carry
51      , RealType* x
52      , RealType modulus)
53  {
54    RealType value;
55    for(unsigned int j = 0; j < f.long_lag; ++j) {
56      is >> value >> std::ws;
57      x[j] = value / modulus;
58    }
59    is >> value >> std::ws;
60    carry = value / modulus;
61  }
62}
63# endif
64// subtract-with-carry generator
65// Marsaglia and Zaman
66
67template<class IntType, IntType m, unsigned int s, unsigned int r,
68  IntType val>
69class subtract_with_carry
70{
71public:
72  typedef IntType result_type;
73  BOOST_STATIC_CONSTANT(bool, has_fixed_range = true);
74  BOOST_STATIC_CONSTANT(result_type, min_value = 0);
75  BOOST_STATIC_CONSTANT(result_type, max_value = m-1);
76  BOOST_STATIC_CONSTANT(result_type, modulus = m);
77  BOOST_STATIC_CONSTANT(unsigned int, long_lag = r);
78  BOOST_STATIC_CONSTANT(unsigned int, short_lag = s);
79
80  subtract_with_carry() {
81    // MSVC fails BOOST_STATIC_ASSERT with std::numeric_limits at class scope
82#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
83    BOOST_STATIC_ASSERT(std::numeric_limits<result_type>::is_signed);
84    BOOST_STATIC_ASSERT(std::numeric_limits<result_type>::is_integer);
85#endif
86    seed();
87  }
88  explicit subtract_with_carry(uint32_t value) { seed(value); }
89  template<class Generator>
90  explicit subtract_with_carry(Generator & gen) { seed(gen); }
91  template<class It> subtract_with_carry(It& first, It last) { seed(first,last); }
92
93  // compiler-generated copy ctor and assignment operator are fine
94
95  void seed(uint32_t value = 19780503u)
96  {
97    random::linear_congruential<int32_t, 40014, 0, 2147483563, 0> intgen(value);
98    seed(intgen);
99  }
100
101  // For GCC, moving this function out-of-line prevents inlining, which may
102  // reduce overall object code size.  However, MSVC does not grok
103  // out-of-line template member functions.
104  template<class Generator>
105  void seed(Generator & gen)
106  {
107    // I could have used std::generate_n, but it takes "gen" by value
108    for(unsigned int j = 0; j < long_lag; ++j)
109      x[j] = gen() % modulus;
110    carry = (x[long_lag-1] == 0);
111    k = 0;
112  }
113
114  template<class It>
115  void seed(It& first, It last)
116  {
117    unsigned int j;
118    for(j = 0; j < long_lag && first != last; ++j, ++first)
119      x[j] = *first % modulus;
120    if(first == last && j < long_lag)
121      throw std::invalid_argument("subtract_with_carry::seed");
122    carry = (x[long_lag-1] == 0);
123    k = 0;
124   }
125
126  result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return min_value; }
127  result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return max_value; }
128
129  result_type operator()()
130  {
131    int short_index = k - short_lag;
132    if(short_index < 0)
133      short_index += long_lag;
134    IntType delta;
135    if (x[short_index] >= x[k] + carry) {
136      // x(n) >= 0
137      delta =  x[short_index] - (x[k] + carry);
138      carry = 0;
139    } else {
140      // x(n) < 0
141      delta = modulus - x[k] - carry + x[short_index];
142      carry = 1;
143    }
144    x[k] = delta;
145    ++k;
146    if(k >= long_lag)
147      k = 0;
148    return delta;
149  }
150
151public:
152  static bool validation(result_type x) { return x == val; }
153 
154#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
155
156#ifndef BOOST_NO_MEMBER_TEMPLATE_FRIENDS
157  template<class CharT, class Traits>
158  friend std::basic_ostream<CharT,Traits>&
159  operator<<(std::basic_ostream<CharT,Traits>& os,
160             const subtract_with_carry& f)
161  {
162    for(unsigned int j = 0; j < f.long_lag; ++j)
163      os << f.compute(j) << " ";
164    os << f.carry << " ";
165    return os;
166  }
167
168  template<class CharT, class Traits>
169  friend std::basic_istream<CharT,Traits>&
170  operator>>(std::basic_istream<CharT,Traits>& is, subtract_with_carry& f)
171  {
172    for(unsigned int j = 0; j < f.long_lag; ++j)
173      is >> f.x[j] >> std::ws;
174    is >> f.carry >> std::ws;
175    f.k = 0;
176    return is;
177  }
178#endif
179
180  friend bool operator==(const subtract_with_carry& x, const subtract_with_carry& y)
181  {
182    for(unsigned int j = 0; j < r; ++j)
183      if(x.compute(j) != y.compute(j))
184        return false;
185    return true;
186  }
187
188  friend bool operator!=(const subtract_with_carry& x, const subtract_with_carry& y)
189  { return !(x == y); }
190#else
191  // Use a member function; Streamable concept not supported.
192  bool operator==(const subtract_with_carry& rhs) const
193  {
194    for(unsigned int j = 0; j < r; ++j)
195      if(compute(j) != rhs.compute(j))
196        return false;
197    return true;
198  }
199
200  bool operator!=(const subtract_with_carry& rhs) const
201  { return !(*this == rhs); }
202#endif
203
204private:
205  // returns x(i-r+index), where index is in 0..r-1
206  IntType compute(unsigned int index) const
207  {
208    return x[(k+index) % long_lag];
209  }
210
211  // state representation; next output (state) is x(i)
212  //   x[0]  ... x[k] x[k+1] ... x[long_lag-1]     represents
213  //  x(i-k) ... x(i) x(i+1) ... x(i-k+long_lag-1)
214  // speed: base: 20-25 nsec
215  // ranlux_4: 230 nsec, ranlux_7: 430 nsec, ranlux_14: 810 nsec
216  // This state representation makes operator== and save/restore more
217  // difficult, because we've already computed "too much" and thus
218  // have to undo some steps to get at x(i-r) etc.
219
220  // state representation: next output (state) is x(i)
221  //   x[0]  ... x[k] x[k+1]          ... x[long_lag-1]     represents
222  //  x(i-k) ... x(i) x(i-long_lag+1) ... x(i-k-1)
223  // speed: base 28 nsec
224  // ranlux_4: 370 nsec, ranlux_7: 688 nsec, ranlux_14: 1343 nsec
225  IntType x[long_lag];
226  unsigned int k;
227  int carry;
228};
229
230#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
231//  A definition is required even for integral static constants
232template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val>
233const bool subtract_with_carry<IntType, m, s, r, val>::has_fixed_range;
234template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val>
235const IntType subtract_with_carry<IntType, m, s, r, val>::min_value;
236template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val>
237const IntType subtract_with_carry<IntType, m, s, r, val>::max_value;
238template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val>
239const IntType subtract_with_carry<IntType, m, s, r, val>::modulus;
240template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val>
241const unsigned int subtract_with_carry<IntType, m, s, r, val>::long_lag;
242template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val>
243const unsigned int subtract_with_carry<IntType, m, s, r, val>::short_lag;
244#endif
245
246
247// use a floating-point representation to produce values in [0..1)
248template<class RealType, int w, unsigned int s, unsigned int r, int val=0>
249class subtract_with_carry_01
250{
251public:
252  typedef RealType result_type;
253  BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
254  BOOST_STATIC_CONSTANT(int, word_size = w);
255  BOOST_STATIC_CONSTANT(unsigned int, long_lag = r);
256  BOOST_STATIC_CONSTANT(unsigned int, short_lag = s);
257
258#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
259  BOOST_STATIC_ASSERT(!std::numeric_limits<result_type>::is_integer);
260#endif
261
262  subtract_with_carry_01() { init_modulus(); seed(); }
263  explicit subtract_with_carry_01(uint32_t value)
264  { init_modulus(); seed(value);   }
265  template<class It> subtract_with_carry_01(It& first, It last)
266  { init_modulus(); seed(first,last); }
267
268private:
269  void init_modulus()
270  {
271#ifndef BOOST_NO_STDC_NAMESPACE
272    // allow for Koenig lookup
273    using std::pow;
274#endif
275    _modulus = pow(RealType(2), word_size);
276  }
277
278public:
279  // compiler-generated copy ctor and assignment operator are fine
280
281  void seed(uint32_t value = 19780503u)
282  {
283#ifndef BOOST_NO_STDC_NAMESPACE
284    // allow for Koenig lookup
285    using std::fmod;
286#endif
287    random::linear_congruential<int32_t, 40014, 0, 2147483563, 0> gen(value);
288    unsigned long array[(w+31)/32 * long_lag];
289    for(unsigned int j = 0; j < sizeof(array)/sizeof(unsigned long); ++j)
290      array[j] = gen();
291    unsigned long * start = array;
292    seed(start, array + sizeof(array)/sizeof(unsigned long));
293  }
294
295  template<class It>
296  void seed(It& first, It last)
297  {
298#ifndef BOOST_NO_STDC_NAMESPACE
299    // allow for Koenig lookup
300    using std::fmod;
301    using std::pow;
302#endif
303    unsigned long mask = ~((~0u) << (w%32));   // now lowest (w%32) bits set
304    RealType two32 = pow(RealType(2), 32);
305    unsigned int j;
306    for(j = 0; j < long_lag && first != last; ++j) {
307      x[j] = RealType(0);
308      for(int i = 0; i < w/32 && first != last; ++i, ++first)
309        x[j] += *first / pow(two32,i+1);
310      if(first != last && mask != 0) {
311        x[j] += fmod((*first & mask) / _modulus, RealType(1));
312        ++first;
313      }
314    }
315    if(first == last && j < long_lag)
316      throw std::invalid_argument("subtract_with_carry_01::seed");
317    carry = (x[long_lag-1] ? 0 : 1 / _modulus);
318    k = 0;
319  }
320
321  result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(0); }
322  result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(1); }
323
324  result_type operator()()
325  {
326    int short_index = k - short_lag;
327    if(short_index < 0)
328      short_index += long_lag;
329    RealType delta = x[short_index] - x[k] - carry;
330    if(delta < 0) {
331      delta += RealType(1);
332      carry = RealType(1)/_modulus;
333    } else {
334      carry = 0;
335    }
336    x[k] = delta;
337    ++k;
338    if(k >= long_lag)
339      k = 0;
340    return delta;
341  }
342
343  static bool validation(result_type x)
344  { return x == val/pow(RealType(2), word_size); }
345 
346#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
347
348#ifndef BOOST_NO_MEMBER_TEMPLATE_FRIENDS
349  template<class CharT, class Traits>
350  friend std::basic_ostream<CharT,Traits>&
351  operator<<(std::basic_ostream<CharT,Traits>& os,
352             const subtract_with_carry_01& f)
353  {
354#ifndef BOOST_NO_STDC_NAMESPACE
355    // allow for Koenig lookup
356    using std::pow;
357#endif
358    std::ios_base::fmtflags oldflags = os.flags(os.dec | os.fixed | os.left);
359    for(unsigned int j = 0; j < f.long_lag; ++j)
360      os << (f.compute(j) * f._modulus) << " ";
361    os << (f.carry * f._modulus);
362    os.flags(oldflags);
363    return os;
364  }
365
366  template<class CharT, class Traits>
367  friend std::basic_istream<CharT,Traits>&
368  operator>>(std::basic_istream<CharT,Traits>& is, subtract_with_carry_01& f)
369  {
370# ifdef BOOST_RANDOM_EXTRACT_SWC_01
371      detail::extract_subtract_with_carry_01(is, f, f.carry, f.x, f._modulus);
372# else
373    // MSVC (up to 7.1) and Borland (up to 5.64) don't handle the template type
374    // parameter "RealType" available from the class template scope, so use
375    // the member typedef
376    typename subtract_with_carry_01::result_type value;
377    for(unsigned int j = 0; j < long_lag; ++j) {
378      is >> value >> std::ws;
379      f.x[j] = value / f._modulus;
380    }
381    is >> value >> std::ws;
382    f.carry = value / f._modulus;
383# endif
384    f.k = 0;
385    return is;
386  }
387#endif
388
389  friend bool operator==(const subtract_with_carry_01& x,
390                         const subtract_with_carry_01& y)
391  {
392    for(unsigned int j = 0; j < r; ++j)
393      if(x.compute(j) != y.compute(j))
394        return false;
395    return true;
396  }
397
398  friend bool operator!=(const subtract_with_carry_01& x,
399                         const subtract_with_carry_01& y)
400  { return !(x == y); }
401#else
402  // Use a member function; Streamable concept not supported.
403  bool operator==(const subtract_with_carry_01& rhs) const
404  {
405    for(unsigned int j = 0; j < r; ++j)
406      if(compute(j) != rhs.compute(j))
407        return false;
408    return true;
409  }
410
411  bool operator!=(const subtract_with_carry_01& rhs) const
412  { return !(*this == rhs); }
413#endif
414
415private:
416  RealType compute(unsigned int index) const;
417  unsigned int k;
418  RealType carry;
419  RealType x[long_lag];
420  RealType _modulus;
421};
422
423#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
424//  A definition is required even for integral static constants
425template<class RealType, int w, unsigned int s, unsigned int r, int val>
426const bool subtract_with_carry_01<RealType, w, s, r, val>::has_fixed_range;
427template<class RealType, int w, unsigned int s, unsigned int r, int val>
428const int subtract_with_carry_01<RealType, w, s, r, val>::word_size;
429template<class RealType, int w, unsigned int s, unsigned int r, int val>
430const unsigned int subtract_with_carry_01<RealType, w, s, r, val>::long_lag;
431template<class RealType, int w, unsigned int s, unsigned int r, int val>
432const unsigned int subtract_with_carry_01<RealType, w, s, r, val>::short_lag;
433#endif
434
435template<class RealType, int w, unsigned int s, unsigned int r, int val>
436RealType subtract_with_carry_01<RealType, w, s, r, val>::compute(unsigned int index) const
437{
438  return x[(k+index) % long_lag];
439}
440
441
442} // namespace random
443} // namespace boost
444
445#endif // BOOST_RANDOM_SUBTRACT_WITH_CARRY_HPP
Note: See TracBrowser for help on using the repository browser.