1 /* boost random/subtract_with_carry.hpp header file
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)
8 * See http://www.boost.org for most recent version including documentation.
10 * $Id: subtract_with_carry.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $
16 #ifndef BOOST_RANDOM_SUBTRACT_WITH_CARRY_HPP
17 #define BOOST_RANDOM_SUBTRACT_WITH_CARRY_HPP
19 #include <boost/config/no_tr1/cmath.hpp>
21 #include <algorithm> // std::equal
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>
37 #if BOOST_WORKAROUND(_MSC_FULL_VER, BOOST_TESTED_AT(13102292)) && BOOST_MSVC > 1300
38 # define BOOST_RANDOM_EXTRACT_SWC_01
41 #if defined(__APPLE_CC__) && defined(__GNUC__) && (__GNUC__ == 3) && (__GNUC_MINOR__ <= 3)
42 # define BOOST_RANDOM_EXTRACT_SWC_01
45 # ifdef BOOST_RANDOM_EXTRACT_SWC_01
48 template <class IStream, class SubtractWithCarry, class RealType>
49 void extract_subtract_with_carry_01(
51 , SubtractWithCarry& f
57 for(unsigned int j = 0; j < f.long_lag; ++j) {
58 is >> value >> std::ws;
59 x[j] = value / modulus;
61 is >> value >> std::ws;
62 carry = value / modulus;
68 * Instantiations of @c subtract_with_carry model a
69 * \pseudo_random_number_generator. The algorithm is
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.
78 template<class IntType, IntType m, unsigned int s, unsigned int r,
80 class subtract_with_carry
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);
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);
99 BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(subtract_with_carry, uint32_t, value)
101 BOOST_RANDOM_DETAIL_GENERATOR_CONSTRUCTOR(subtract_with_carry, Generator, gen)
103 template<class It> subtract_with_carry(It& first, It last) { seed(first,last); }
105 // compiler-generated copy ctor and assignment operator are fine
107 void seed() { seed(19780503u); }
108 BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(subtract_with_carry, uint32_t, value)
110 random::linear_congruential<int32_t, 40014, 0, 2147483563, 0> intgen(value);
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)
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);
127 void seed(It& first, It last)
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);
138 result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return min_value; }
139 result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return max_value; }
141 result_type operator()()
143 int short_index = k - short_lag;
145 short_index += long_lag;
147 if (x[short_index] >= x[k] + carry) {
149 delta = x[short_index] - (x[k] + carry);
153 delta = modulus - x[k] - carry + x[short_index];
164 static bool validation(result_type x) { return x == val; }
166 #ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
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)
174 for(unsigned int j = 0; j < f.long_lag; ++j)
175 os << f.compute(j) << " ";
176 os << f.carry << " ";
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)
184 for(unsigned int j = 0; j < f.long_lag; ++j)
185 is >> f.x[j] >> std::ws;
186 is >> f.carry >> std::ws;
192 friend bool operator==(const subtract_with_carry& x, const subtract_with_carry& y)
194 for(unsigned int j = 0; j < r; ++j)
195 if(x.compute(j) != y.compute(j))
200 friend bool operator!=(const subtract_with_carry& x, const subtract_with_carry& y)
201 { return !(x == y); }
203 // Use a member function; Streamable concept not supported.
204 bool operator==(const subtract_with_carry& rhs) const
206 for(unsigned int j = 0; j < r; ++j)
207 if(compute(j) != rhs.compute(j))
212 bool operator!=(const subtract_with_carry& rhs) const
213 { return !(*this == rhs); }
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
221 return x[(k+index) % long_lag];
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.
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
244 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
245 // A definition is required even for integral static constants
246 template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val>
247 const bool subtract_with_carry<IntType, m, s, r, val>::has_fixed_range;
248 template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val>
249 const IntType subtract_with_carry<IntType, m, s, r, val>::min_value;
250 template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val>
251 const IntType subtract_with_carry<IntType, m, s, r, val>::max_value;
252 template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val>
253 const IntType subtract_with_carry<IntType, m, s, r, val>::modulus;
254 template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val>
255 const unsigned int subtract_with_carry<IntType, m, s, r, val>::long_lag;
256 template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val>
257 const unsigned int subtract_with_carry<IntType, m, s, r, val>::short_lag;
261 // use a floating-point representation to produce values in [0..1)
262 /** @copydoc boost::random::subtract_with_carry */
263 template<class RealType, int w, unsigned int s, unsigned int r, int val=0>
264 class subtract_with_carry_01
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);
273 #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
274 BOOST_STATIC_ASSERT(!std::numeric_limits<result_type>::is_integer);
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); }
284 /// \cond hide_private_members
287 #ifndef BOOST_NO_STDC_NAMESPACE
288 // allow for Koenig lookup
291 _modulus = pow(RealType(2), word_size);
293 /// \endcond hide_private_members
296 // compiler-generated copy ctor and assignment operator are fine
298 void seed(uint32_t value = 19780503u)
300 #ifndef BOOST_NO_STDC_NAMESPACE
301 // allow for Koenig lookup
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)
308 unsigned long * start = array;
309 seed(start, array + sizeof(array)/sizeof(unsigned long));
313 void seed(It& first, It last)
315 #ifndef BOOST_NO_STDC_NAMESPACE
316 // allow for Koenig lookup
320 unsigned long mask = ~((~0u) << (w%32)); // now lowest (w%32) bits set
321 RealType two32 = pow(RealType(2), 32);
323 for(j = 0; j < long_lag && first != last; ++j) {
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));
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);
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); }
341 result_type operator()()
343 int short_index = k - short_lag;
345 short_index += long_lag;
346 RealType delta = x[short_index] - x[k] - carry;
348 delta += RealType(1);
349 carry = RealType(1)/_modulus;
360 static bool validation(result_type x)
361 { return x == val/pow(RealType(2), word_size); }
363 #ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
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)
371 #ifndef BOOST_NO_STDC_NAMESPACE
372 // allow for Koenig lookup
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);
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)
387 # ifdef BOOST_RANDOM_EXTRACT_SWC_01
388 detail::extract_subtract_with_carry_01(is, f, f.carry, f.x, f._modulus);
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;
398 is >> value >> std::ws;
399 f.carry = value / f._modulus;
406 friend bool operator==(const subtract_with_carry_01& x,
407 const subtract_with_carry_01& y)
409 for(unsigned int j = 0; j < r; ++j)
410 if(x.compute(j) != y.compute(j))
415 friend bool operator!=(const subtract_with_carry_01& x,
416 const subtract_with_carry_01& y)
417 { return !(x == y); }
419 // Use a member function; Streamable concept not supported.
420 bool operator==(const subtract_with_carry_01& rhs) const
422 for(unsigned int j = 0; j < r; ++j)
423 if(compute(j) != rhs.compute(j))
428 bool operator!=(const subtract_with_carry_01& rhs) const
429 { return !(*this == rhs); }
433 /// \cond hide_private_members
434 RealType compute(unsigned int index) const;
438 RealType x[long_lag];
442 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
443 // A definition is required even for integral static constants
444 template<class RealType, int w, unsigned int s, unsigned int r, int val>
445 const bool subtract_with_carry_01<RealType, w, s, r, val>::has_fixed_range;
446 template<class RealType, int w, unsigned int s, unsigned int r, int val>
447 const int subtract_with_carry_01<RealType, w, s, r, val>::word_size;
448 template<class RealType, int w, unsigned int s, unsigned int r, int val>
449 const unsigned int subtract_with_carry_01<RealType, w, s, r, val>::long_lag;
450 template<class RealType, int w, unsigned int s, unsigned int r, int val>
451 const unsigned int subtract_with_carry_01<RealType, w, s, r, val>::short_lag;
454 /// \cond hide_private_members
455 template<class RealType, int w, unsigned int s, unsigned int r, int val>
456 RealType subtract_with_carry_01<RealType, w, s, r, val>::compute(unsigned int index) const
458 return x[(k+index) % long_lag];
462 } // namespace random
465 #endif // BOOST_RANDOM_SUBTRACT_WITH_CARRY_HPP