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$
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> // std::pow
20#include <iostream>
21#include <algorithm> // std::equal
22#include <stdexcept>
23#include <boost/config.hpp>
24#include <boost/limits.hpp>
25#include <boost/cstdint.hpp>
26#include <boost/static_assert.hpp>
27#include <boost/integer/static_log2.hpp>
28#include <boost/integer/integer_mask.hpp>
29#include <boost/detail/workaround.hpp>
30#include <boost/random/detail/config.hpp>
31#include <boost/random/detail/seed.hpp>
32#include <boost/random/detail/operators.hpp>
33#include <boost/random/detail/seed_impl.hpp>
34#include <boost/random/detail/generator_seed_seq.hpp>
35#include <boost/random/linear_congruential.hpp>
36
37
38namespace boost {
39namespace random {
40
41namespace detail {
42
43struct subtract_with_carry_discard
44{
45 template<class Engine>
46 static void apply(Engine& eng, boost::uintmax_t z)
47 {
48 typedef typename Engine::result_type IntType;
49 const std::size_t short_lag = Engine::short_lag;
50 const std::size_t long_lag = Engine::long_lag;
51 std::size_t k = eng.k;
52 IntType carry = eng.carry;
53 if(k != 0) {
54 // increment k until it becomes 0.
55 if(k < short_lag) {
56 std::size_t limit = (short_lag - k) < z?
57 short_lag : (k + static_cast<std::size_t>(z));
58 for(std::size_t j = k; j < limit; ++j) {
59 carry = eng.do_update(j, j + long_lag - short_lag, carry);
60 }
61 }
62 std::size_t limit = (long_lag - k) < z?
63 long_lag : (k + static_cast<std::size_t>(z));
64 std::size_t start = (k < short_lag ? short_lag : k);
65 for(std::size_t j = start; j < limit; ++j) {
66 carry = eng.do_update(j, j - short_lag, carry);
67 }
68 }
69
70 k = ((z % long_lag) + k) % long_lag;
71
72 if(k < z) {
73 // main loop: update full blocks from k = 0 to long_lag
74 for(std::size_t i = 0; i < (z - k) / long_lag; ++i) {
75 for(std::size_t j = 0; j < short_lag; ++j) {
76 carry = eng.do_update(j, j + long_lag - short_lag, carry);
77 }
78 for(std::size_t j = short_lag; j < long_lag; ++j) {
79 carry = eng.do_update(j, j - short_lag, carry);
80 }
81 }
82
83 // Update the last partial block
84 std::size_t limit = short_lag < k? short_lag : k;
85 for(std::size_t j = 0; j < limit; ++j) {
86 carry = eng.do_update(j, j + long_lag - short_lag, carry);
87 }
88 for(std::size_t j = short_lag; j < k; ++j) {
89 carry = eng.do_update(j, j - short_lag, carry);
90 }
91 }
92 eng.carry = carry;
93 eng.k = k;
94 }
95};
96
97}
98
99/**
100 * Instantiations of @c subtract_with_carry_engine model a
101 * \pseudo_random_number_generator. The algorithm is
102 * described in
103 *
104 * @blockquote
105 * "A New Class of Random Number Generators", George
106 * Marsaglia and Arif Zaman, Annals of Applied Probability,
107 * Volume 1, Number 3 (1991), 462-480.
108 * @endblockquote
109 */
110template<class IntType, std::size_t w, std::size_t s, std::size_t r>
111class subtract_with_carry_engine
112{
113public:
114 typedef IntType result_type;
115 BOOST_STATIC_CONSTANT(std::size_t, word_size = w);
116 BOOST_STATIC_CONSTANT(std::size_t, long_lag = r);
117 BOOST_STATIC_CONSTANT(std::size_t, short_lag = s);
118 BOOST_STATIC_CONSTANT(uint32_t, default_seed = 19780503u);
119
120 // Required by the old Boost.Random concepts
121 BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
122 // Backwards compatibility
123 BOOST_STATIC_CONSTANT(result_type, modulus = (result_type(1) << w));
124
125 BOOST_STATIC_ASSERT(std::numeric_limits<result_type>::is_integer);
126
127 /**
128 * Constructs a new @c subtract_with_carry_engine and seeds
129 * it with the default seed.
130 */
131 subtract_with_carry_engine() { seed(); }
132 /**
133 * Constructs a new @c subtract_with_carry_engine and seeds
134 * it with @c value.
135 */
136 BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(subtract_with_carry_engine,
137 IntType, value)
138 { seed(value); }
139 /**
140 * Constructs a new @c subtract_with_carry_engine and seeds
141 * it with values produced by @c seq.generate().
142 */
143 BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(subtract_with_carry_engine,
144 SeedSeq, seq)
145 { seed(seq); }
146 /**
147 * Constructs a new @c subtract_with_carry_engine and seeds
148 * it with values from a range. first is updated to point
149 * one past the last value consumed. If there are not
150 * enough elements in the range to fill the entire state of
151 * the generator, throws @c std::invalid_argument.
152 */
153 template<class It> subtract_with_carry_engine(It& first, It last)
154 { seed(first,last); }
155
156 // compiler-generated copy ctor and assignment operator are fine
157
158 /** Seeds the generator with the default seed. */
159 void seed() { seed(default_seed); }
160 BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(subtract_with_carry_engine,
161 IntType, value)
162 {
163 typedef linear_congruential_engine<uint32_t,40014,0,2147483563> gen_t;
164 gen_t intgen(static_cast<boost::uint32_t>(value == 0 ? default_seed : value));
165 detail::generator_seed_seq<gen_t> gen(intgen);
166 seed(gen);
167 }
168
169 /** Seeds the generator with values produced by @c seq.generate(). */
170 BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(subtract_with_carry, SeedSeq, seq)
171 {
172 detail::seed_array_int<w>(seq, x);
173 carry = (x[long_lag-1] == 0);
174 k = 0;
175 }
176
177 /**
178 * Seeds the generator with values from a range. Updates @c first to
179 * point one past the last consumed value. If the range does not
180 * contain enough elements to fill the entire state of the generator,
181 * throws @c std::invalid_argument.
182 */
183 template<class It>
184 void seed(It& first, It last)
185 {
186 detail::fill_array_int<w>(first, last, x);
187 carry = (x[long_lag-1] == 0);
188 k = 0;
189 }
190
191 /** Returns the smallest value that the generator can produce. */
192 static result_type min BOOST_PREVENT_MACRO_SUBSTITUTION ()
193 { return 0; }
194 /** Returns the largest value that the generator can produce. */
195 static result_type max BOOST_PREVENT_MACRO_SUBSTITUTION ()
196 { return boost::low_bits_mask_t<w>::sig_bits; }
197
198 /** Returns the next value of the generator. */
199 result_type operator()()
200 {
201 std::size_t short_index =
202 (k < short_lag)?
203 (k + long_lag - short_lag) :
204 (k - short_lag);
205 carry = do_update(k, short_index, carry);
206 IntType result = x[k];
207 ++k;
208 if(k >= long_lag)
209 k = 0;
210 return result;
211 }
212
213 /** Advances the state of the generator by @c z. */
214 void discard(boost::uintmax_t z)
215 {
216 detail::subtract_with_carry_discard::apply(*this, z);
217 }
218
219 /** Fills a range with random values. */
220 template<class It>
221 void generate(It first, It last)
222 { detail::generate_from_int(*this, first, last); }
223
224 /** Writes a @c subtract_with_carry_engine to a @c std::ostream. */
225 BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, subtract_with_carry_engine, f)
226 {
227 for(unsigned int j = 0; j < f.long_lag; ++j)
228 os << f.compute(j) << ' ';
229 os << f.carry;
230 return os;
231 }
232
233 /** Reads a @c subtract_with_carry_engine from a @c std::istream. */
234 BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, subtract_with_carry_engine, f)
235 {
236 for(unsigned int j = 0; j < f.long_lag; ++j)
237 is >> f.x[j] >> std::ws;
238 is >> f.carry;
239 f.k = 0;
240 return is;
241 }
242
243 /**
244 * Returns true if the two generators will produce identical
245 * sequences of values.
246 */
247 BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(subtract_with_carry_engine, x, y)
248 {
249 for(unsigned int j = 0; j < r; ++j)
250 if(x.compute(j) != y.compute(j))
251 return false;
252 return true;
253 }
254
255 /**
256 * Returns true if the two generators will produce different
257 * sequences of values.
258 */
259 BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(subtract_with_carry_engine)
260
261private:
262 /// \cond show_private
263 // returns x(i-r+index), where index is in 0..r-1
264 IntType compute(unsigned int index) const
265 {
266 return x[(k+index) % long_lag];
267 }
268
269 friend struct detail::subtract_with_carry_discard;
270
271 IntType do_update(std::size_t current, std::size_t short_index, IntType carry)
272 {
273 IntType delta;
274 IntType temp = x[current] + carry;
275 if (x[short_index] >= temp) {
276 // x(n) >= 0
277 delta = x[short_index] - temp;
278 carry = 0;
279 } else {
280 // x(n) < 0
281 delta = modulus - temp + x[short_index];
282 carry = 1;
283 }
284 x[current] = delta;
285 return carry;
286 }
287 /// \endcond
288
289 // state representation; next output (state) is x(i)
290 // x[0] ... x[k] x[k+1] ... x[long_lag-1] represents
291 // x(i-k) ... x(i) x(i+1) ... x(i-k+long_lag-1)
292 // speed: base: 20-25 nsec
293 // ranlux_4: 230 nsec, ranlux_7: 430 nsec, ranlux_14: 810 nsec
294 // This state representation makes operator== and save/restore more
295 // difficult, because we've already computed "too much" and thus
296 // have to undo some steps to get at x(i-r) etc.
297
298 // state representation: next output (state) is x(i)
299 // x[0] ... x[k] x[k+1] ... x[long_lag-1] represents
300 // x(i-k) ... x(i) x(i-long_lag+1) ... x(i-k-1)
301 // speed: base 28 nsec
302 // ranlux_4: 370 nsec, ranlux_7: 688 nsec, ranlux_14: 1343 nsec
303 IntType x[long_lag];
304 std::size_t k;
305 IntType carry;
306};
307
308#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
309// A definition is required even for integral static constants
310template<class IntType, std::size_t w, std::size_t s, std::size_t r>
311const bool subtract_with_carry_engine<IntType, w, s, r>::has_fixed_range;
312template<class IntType, std::size_t w, std::size_t s, std::size_t r>
313const IntType subtract_with_carry_engine<IntType, w, s, r>::modulus;
314template<class IntType, std::size_t w, std::size_t s, std::size_t r>
315const std::size_t subtract_with_carry_engine<IntType, w, s, r>::word_size;
316template<class IntType, std::size_t w, std::size_t s, std::size_t r>
317const std::size_t subtract_with_carry_engine<IntType, w, s, r>::long_lag;
318template<class IntType, std::size_t w, std::size_t s, std::size_t r>
319const std::size_t subtract_with_carry_engine<IntType, w, s, r>::short_lag;
320template<class IntType, std::size_t w, std::size_t s, std::size_t r>
321const uint32_t subtract_with_carry_engine<IntType, w, s, r>::default_seed;
322#endif
323
324
325// use a floating-point representation to produce values in [0..1)
326/**
327 * Instantiations of \subtract_with_carry_01_engine model a
328 * \pseudo_random_number_generator. The algorithm is
329 * described in
330 *
331 * @blockquote
332 * "A New Class of Random Number Generators", George
333 * Marsaglia and Arif Zaman, Annals of Applied Probability,
334 * Volume 1, Number 3 (1991), 462-480.
335 * @endblockquote
336 */
337template<class RealType, std::size_t w, std::size_t s, std::size_t r>
338class subtract_with_carry_01_engine
339{
340public:
341 typedef RealType result_type;
342 BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
343 BOOST_STATIC_CONSTANT(std::size_t, word_size = w);
344 BOOST_STATIC_CONSTANT(std::size_t, long_lag = r);
345 BOOST_STATIC_CONSTANT(std::size_t, short_lag = s);
346 BOOST_STATIC_CONSTANT(boost::uint32_t, default_seed = 19780503u);
347
348 BOOST_STATIC_ASSERT(!std::numeric_limits<result_type>::is_integer);
349
350 /** Creates a new \subtract_with_carry_01_engine using the default seed. */
351 subtract_with_carry_01_engine() { init_modulus(); seed(); }
352 /** Creates a new subtract_with_carry_01_engine and seeds it with value. */
353 BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(subtract_with_carry_01_engine,
354 boost::uint32_t, value)
355 { init_modulus(); seed(value); }
356 /**
357 * Creates a new \subtract_with_carry_01_engine and seeds with values
358 * produced by seq.generate().
359 */
360 BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(subtract_with_carry_01_engine,
361 SeedSeq, seq)
362 { init_modulus(); seed(seq); }
363 /**
364 * Creates a new \subtract_with_carry_01_engine and seeds it with values
365 * from a range. Advances first to point one past the last consumed
366 * value. If the range does not contain enough elements to fill the
367 * entire state, throws @c std::invalid_argument.
368 */
369 template<class It> subtract_with_carry_01_engine(It& first, It last)
370 { init_modulus(); seed(first,last); }
371
372private:
373 /// \cond show_private
374 void init_modulus()
375 {
376#ifndef BOOST_NO_STDC_NAMESPACE
377 // allow for Koenig lookup
378 using std::pow;
379#endif
380 _modulus = pow(RealType(2), RealType(word_size));
381 }
382 /// \endcond
383
384public:
385 // compiler-generated copy ctor and assignment operator are fine
386
387 /** Seeds the generator with the default seed. */
388 void seed() { seed(default_seed); }
389
390 /** Seeds the generator with @c value. */
391 BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(subtract_with_carry_01_engine,
392 boost::uint32_t, value)
393 {
394 typedef linear_congruential_engine<uint32_t, 40014, 0, 2147483563> gen_t;
395 gen_t intgen(value == 0 ? default_seed : value);
396 detail::generator_seed_seq<gen_t> gen(intgen);
397 seed(gen);
398 }
399
400 /** Seeds the generator with values produced by @c seq.generate(). */
401 BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(subtract_with_carry_01_engine,
402 SeedSeq, seq)
403 {
404 detail::seed_array_real<w>(seq, x);
405 carry = (x[long_lag-1] ? result_type(0) : result_type(1 / _modulus));
406 k = 0;
407 }
408
409 /**
410 * Seeds the generator with values from a range. Updates first to
411 * point one past the last consumed element. If there are not
412 * enough elements in the range to fill the entire state, throws
413 * @c std::invalid_argument.
414 */
415 template<class It>
416 void seed(It& first, It last)
417 {
418 detail::fill_array_real<w>(first, last, x);
419 carry = (x[long_lag-1] ? result_type(0) : result_type(1 / _modulus));
420 k = 0;
421 }
422
423 /** Returns the smallest value that the generator can produce. */
424 static result_type min BOOST_PREVENT_MACRO_SUBSTITUTION ()
425 { return result_type(0); }
426 /** Returns the largest value that the generator can produce. */
427 static result_type max BOOST_PREVENT_MACRO_SUBSTITUTION ()
428 { return result_type(1); }
429
430 /** Returns the next value of the generator. */
431 result_type operator()()
432 {
433 std::size_t short_index =
434 (k < short_lag) ?
435 (k + long_lag - short_lag) :
436 (k - short_lag);
437 carry = do_update(k, short_index, carry);
438 RealType result = x[k];
439 ++k;
440 if(k >= long_lag)
441 k = 0;
442 return result;
443 }
444
445 /** Advances the state of the generator by @c z. */
446 void discard(boost::uintmax_t z)
447 { detail::subtract_with_carry_discard::apply(*this, z); }
448
449 /** Fills a range with random values. */
450 template<class Iter>
451 void generate(Iter first, Iter last)
452 { detail::generate_from_real(*this, first, last); }
453
454 /** Writes a \subtract_with_carry_01_engine to a @c std::ostream. */
455 BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, subtract_with_carry_01_engine, f)
456 {
457 std::ios_base::fmtflags oldflags =
458 os.flags(os.dec | os.fixed | os.left);
459 for(unsigned int j = 0; j < f.long_lag; ++j)
460 os << (f.compute(j) * f._modulus) << ' ';
461 os << (f.carry * f._modulus);
462 os.flags(oldflags);
463 return os;
464 }
465
466 /** Reads a \subtract_with_carry_01_engine from a @c std::istream. */
467 BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, subtract_with_carry_01_engine, f)
468 {
469 RealType value;
470 for(unsigned int j = 0; j < long_lag; ++j) {
471 is >> value >> std::ws;
472 f.x[j] = value / f._modulus;
473 }
474 is >> value;
475 f.carry = value / f._modulus;
476 f.k = 0;
477 return is;
478 }
479
480 /** Returns true if the two generators will produce identical sequences. */
481 BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(subtract_with_carry_01_engine, x, y)
482 {
483 for(unsigned int j = 0; j < r; ++j)
484 if(x.compute(j) != y.compute(j))
485 return false;
486 return true;
487 }
488
489 /** Returns true if the two generators will produce different sequences. */
490 BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(subtract_with_carry_01_engine)
491
492private:
493 /// \cond show_private
494 RealType compute(unsigned int index) const
495 {
496 return x[(k+index) % long_lag];
497 }
498
499 friend struct detail::subtract_with_carry_discard;
500
501 RealType do_update(std::size_t current, std::size_t short_index, RealType carry)
502 {
503 RealType delta = x[short_index] - x[current] - carry;
504 if(delta < 0) {
505 delta += RealType(1);
506 carry = RealType(1)/_modulus;
507 } else {
508 carry = 0;
509 }
510 x[current] = delta;
511 return carry;
512 }
513 /// \endcond
514 std::size_t k;
515 RealType carry;
516 RealType x[long_lag];
517 RealType _modulus;
518};
519
520#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
521// A definition is required even for integral static constants
522template<class RealType, std::size_t w, std::size_t s, std::size_t r>
523const bool subtract_with_carry_01_engine<RealType, w, s, r>::has_fixed_range;
524template<class RealType, std::size_t w, std::size_t s, std::size_t r>
525const std::size_t subtract_with_carry_01_engine<RealType, w, s, r>::word_size;
526template<class RealType, std::size_t w, std::size_t s, std::size_t r>
527const std::size_t subtract_with_carry_01_engine<RealType, w, s, r>::long_lag;
528template<class RealType, std::size_t w, std::size_t s, std::size_t r>
529const std::size_t subtract_with_carry_01_engine<RealType, w, s, r>::short_lag;
530template<class RealType, std::size_t w, std::size_t s, std::size_t r>
531const uint32_t subtract_with_carry_01_engine<RealType, w, s, r>::default_seed;
532#endif
533
534
535/// \cond show_deprecated
536
537template<class IntType, IntType m, unsigned s, unsigned r, IntType v>
538class subtract_with_carry :
539 public subtract_with_carry_engine<IntType,
540 boost::static_log2<m>::value, s, r>
541{
542 typedef subtract_with_carry_engine<IntType,
543 boost::static_log2<m>::value, s, r> base_type;
544public:
545 subtract_with_carry() {}
546 BOOST_RANDOM_DETAIL_GENERATOR_CONSTRUCTOR(subtract_with_carry, Gen, gen)
547 { seed(gen); }
548 BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(subtract_with_carry,
549 IntType, val)
550 { seed(val); }
551 template<class It>
552 subtract_with_carry(It& first, It last) : base_type(first, last) {}
553 void seed() { base_type::seed(); }
554 BOOST_RANDOM_DETAIL_GENERATOR_SEED(subtract_with_carry, Gen, gen)
555 {
556 detail::generator_seed_seq<Gen> seq(gen);
557 base_type::seed(seq);
558 }
559 BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(subtract_with_carry, IntType, val)
560 { base_type::seed(val); }
561 template<class It>
562 void seed(It& first, It last) { base_type::seed(first, last); }
563};
564
565template<class RealType, int w, unsigned s, unsigned r, int v = 0>
566class subtract_with_carry_01 :
567 public subtract_with_carry_01_engine<RealType, w, s, r>
568{
569 typedef subtract_with_carry_01_engine<RealType, w, s, r> base_type;
570public:
571 subtract_with_carry_01() {}
572 BOOST_RANDOM_DETAIL_GENERATOR_CONSTRUCTOR(subtract_with_carry_01, Gen, gen)
573 { seed(gen); }
574 BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(subtract_with_carry_01,
575 uint32_t, val)
576 { seed(val); }
577 template<class It>
578 subtract_with_carry_01(It& first, It last) : base_type(first, last) {}
579 void seed() { base_type::seed(); }
580 BOOST_RANDOM_DETAIL_GENERATOR_SEED(subtract_with_carry_01, Gen, gen)
581 {
582 detail::generator_seed_seq<Gen> seq(gen);
583 base_type::seed(seq);
584 }
585 BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(subtract_with_carry_01, uint32_t, val)
586 { base_type::seed(val); }
587 template<class It>
588 void seed(It& first, It last) { base_type::seed(first, last); }
589};
590
591/// \endcond
592
593namespace detail {
594
595template<class Engine>
596struct generator_bits;
597
598template<class RealType, std::size_t w, std::size_t s, std::size_t r>
599struct generator_bits<subtract_with_carry_01_engine<RealType, w, s, r> > {
600 static std::size_t value() { return w; }
601};
602
603template<class RealType, int w, unsigned s, unsigned r, int v>
604struct generator_bits<subtract_with_carry_01<RealType, w, s, r, v> > {
605 static std::size_t value() { return w; }
606};
607
608}
609
610} // namespace random
611} // namespace boost
612
613#endif // BOOST_RANDOM_SUBTRACT_WITH_CARRY_HPP
614