PageRenderTime 56ms CodeModel.GetById 20ms app.highlight 30ms RepoModel.GetById 1ms app.codeStats 0ms

/src/contrib/boost/random/subtract_with_carry.hpp

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