RandomLib  1.10
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
RandomNumber.hpp
Go to the documentation of this file.
1 /**
2  * \file RandomNumber.hpp
3  * \brief Header for RandomNumber
4  *
5  * Infinite precision random numbers.
6  *
7  * Copyright (c) Charles Karney (2006-2013) <charles@karney.com> and licensed
8  * under the MIT/X11 License. For more information, see
9  * http://randomlib.sourceforge.net/
10  **********************************************************************/
11 
12 #if !defined(RANDOMLIB_RANDOMNUMBER_HPP)
13 #define RANDOMLIB_RANDOMNUMBER_HPP 1
14 
15 #include <vector>
16 #include <iomanip>
17 #include <limits>
18 #include <cmath> // for std::pow
20 
21 namespace RandomLib {
22  /**
23  * \brief Infinite precision random numbers.
24  *
25  * Implement infinite precision random numbers. Integer part is non-random.
26  * Fraction part consists of any some number of digits in base
27  * 2<sup><i>b</i></sup>. If \e m digits have been generated then the
28  * fraction is uniformly distributed in the open interval
29  * &sum;<sub><i>k</i>=1</sub><sup><i>m</i></sup>
30  * <i>f</i><sub><i>k</i>&minus;1</sub>/2<sup><i>kb</i></sup> +
31  * (0,1)/2<sup><i>mb</i></sup>. When a RandomNumber is first constructed the
32  * integer part is zero and \e m = 0, and the number represents (0,1). A
33  * RandomNumber is able to represent all numbers in the symmetric open
34  * interval (&minus;2<sup>31</sup>, 2<sup>31</sup>). In this implementation,
35  * \e b must one of 1, 2, 3, 4, 8, 12, 16, 20, 24, 28, or 32. (This
36  * restriction allows printing in hexadecimal and can easily be relaxed.
37  * There's also no essential reason why the base should be a power of 2.)
38  *
39  * @tparam bits the number of bits in each digit.
40  **********************************************************************/
41  template<int bits = 1> class RandomNumber {
42  public:
43  /**
44  * Constructor sets number to a random number uniformly distributed in
45  * (0,1).
46  **********************************************************************/
47  RandomNumber() throw() : _n(0), _s(1) {}
48  /**
49  * Swap with another RandomNumber. This is a fast way of doing an
50  * assignment.
51  *
52  * @param[in,out] t the RandomNumber to swap with.
53  **********************************************************************/
54  void swap(RandomNumber& t) throw() {
55  if (this != &t) {
56  std::swap(_n, t._n);
57  std::swap(_s, t._s);
58  _f.swap(t._f);
59  }
60  }
61  /**
62  * Return to initial state, uniformly distributed in (0,1).
63  **********************************************************************/
64  void Init() throw() {
65  STATIC_ASSERT(bits > 0 && bits <= w && (bits < 4 || bits % 4 == 0),
66  "RandomNumber: unsupported value for bits");
67  _n = 0;
68  _s = 1;
69  _f.clear();
70  }
71  /**
72  * @return the sign of the RandomNumber (&plusmn; 1).
73  **********************************************************************/
74  int Sign() const throw() { return _s; }
75  /**
76  * Change the sign of the RandomNumber.
77  **********************************************************************/
78  void Negate() throw() { _s *= -1; }
79  /**
80  * @return the floor of the RandomNumber.
81  **********************************************************************/
82  int Floor() const throw() { return _s > 0 ? int(_n) : -1 - int(_n); }
83  /**
84  * @return the ceiling of the RandomNumber.
85  **********************************************************************/
86  int Ceiling() const throw() { return _s > 0 ? 1 + int(_n) : - int(_n); }
87  /**
88  * @return the unsigned integer component of the RandomNumber.
89  **********************************************************************/
90  unsigned UInteger() const throw() { return _n; }
91  /**
92  * Add integer \e k to the RandomNumber.
93  *
94  * @param[in] k the integer to add.
95  **********************************************************************/
96  void AddInteger(int k) throw() {
97  k += Floor(); // The new floor
98  int ns = k < 0 ? -1 : 1; // The new sign
99  if (ns != _s) // If sign changes, set f = 1 - f
100  for (size_t k = 0; k < Size(); ++k)
101  _f[k] = ~_f[k] & mask;
102  _n = ns > 0 ? k : -(k + 1);
103  }
104  /**
105  * Compare with another RandomNumber, *this &lt; \e t
106  *
107  * @tparam Random the type of the random generator.
108  * @param[in,out] r a random generator.
109  * @param[in,out] t a RandomNumber to compare.
110  * @return true if *this &lt; \e t.
111  **********************************************************************/
112  template<class Random> bool LessThan(Random& r, RandomNumber& t) {
113  if (this == &t) return false; // same object
114  if (_s != t._s) return _s < t._s;
115  if (_n != t._n) return (_s < 0) ^ (_n < t._n);
116  for (unsigned k = 0; ; ++k) {
117  // Impose an order on the evaluation of the digits.
118  const unsigned x = Digit(r,k);
119  const unsigned y = t.Digit(r,k);
120  if (x != y) return (_s < 0) ^ (x < y);
121  // Two distinct numbers are never equal
122  }
123  }
124  /**
125  * Compare RandomNumber with two others, *this &gt; max(\e u, \e v)
126  *
127  * @tparam Random the type of the random generator.
128  * @param[in,out] r a random generator.
129  * @param[in,out] u first RandomNumber to compare.
130  * @param[in,out] v second RandomNumber to compare.
131  * @return true if *this &gt; max(\e u, \e v).
132  **********************************************************************/
133  template<class Random> bool GreaterPair(Random& r,
134  RandomNumber& u, RandomNumber& v) {
135  // cmps is set to false as soon as u <= *this, and likewise for cmpt.
136  bool cmpu = this != &u, cmpv = this != &v && &u != &v;
137  if (!(cmpu || cmpv)) return true;
138  // Check signs first
139  if (cmpu) {
140  if (u._s > _s) return false; // u > *this
141  if (u._s < _s) cmpu = false;
142  }
143  if (cmpv) {
144  if (v._s > _s) return false; // v > *this
145  if (v._s < _s) cmpv = false;
146  }
147  if (!(cmpu || cmpv)) return true; // u <= *this && v <= *this
148  // Check integer parts
149  if (cmpu) {
150  if ((_s < 0) ^ (u._n > _n)) return false; // u > *this
151  if ((_s < 0) ^ (u._n < _n)) cmpu = false;
152  }
153  if (cmpv) {
154  if ((_s < 0) ^ (v._n > _n)) return false; // v > *this
155  if ((_s < 0) ^ (v._n < _n)) cmpv = false;
156  }
157  if (!(cmpu || cmpv)) return true; // u <= *this && v <= *this
158  // Check fractions
159  for (unsigned k = 0; ; ++k) {
160  // Impose an order on the evaluation of the digits. Note that this is
161  // asymmetric on interchange of u and v; since u is tested first, more
162  // digits of u are generated than v (on average).
163  const unsigned x = Digit(r,k);
164  if (cmpu) {
165  const unsigned y = u.Digit(r,k);
166  if ((_s < 0) ^ (y > x)) return false; // u > *this
167  if ((_s < 0) ^ (y < x)) cmpu = false;
168  }
169  if (cmpv) {
170  const unsigned y = v.Digit(r,k);
171  if ((_s < 0) ^ (y > x)) return false; // v > *this
172  if ((_s < 0) ^ (y < x)) cmpv = false;
173  }
174  if (!(cmpu || cmpv)) return true; // u <= *this && v <= *this
175  }
176  }
177  /**
178  * Compare with a fraction, *this &lt; <i>p</i>/<i>q</i>
179  *
180  * @tparam Random the type of the random generator.
181  * @param[in,out] r a random generator.
182  * @param[in] p the numerator of the fraction.
183  * @param[in] q the denominator of the fraction (require \e q &gt; 0).
184  * @return true if *this &lt; <i>p</i>/<i>q</i>.
185  **********************************************************************/
186  template<class Random, typename IntType>
187  bool LessThan(Random& r, IntType p, IntType q) {
188  for (int k = 0;; ++k) {
189  if (p <= 0) return false;
190  if (p >= q) return true;
191  // Here p is in [1,q-1]. Need to avoid overflow in computation of
192  // (q-1)<<bits and (2^bits-1)*q
193  p = (p << bits) - Digit(r,k) * q;
194  }
195  }
196  /**
197  * Compare with a paritally sampled fraction
198  *
199  * @tparam Random the type of the random generator.
200  * @param[in,out] r a random generator.
201  * @param[in] p0 the starting point for the numerator.
202  * @param[in] c the stride for the fraction (require \e c &gt; 0).
203  * @param[in] q the denominator of the fraction (require \e q &gt; 0).
204  * @param[in,out] j the increment for the numerator.
205  * @return true if *this &lt; (<i>p</i><sub>0</sub> + <i>cj</i>)/<i>q</i>.
206  **********************************************************************/
207  template<class Random, typename IntType>
208  bool LessThan(Random& r, IntType p0, IntType c, IntType q,
210  for (int k = 0;; ++k) {
211  if (j. LessThanEqual(r, - p0, c)) return false;
212  if (j.GreaterThanEqual(r, q - p0, c)) return true;
213  p0 = (p0 << bits) - IntType(Digit(r,k)) * q;
214  c <<= bits;
215  }
216  }
217 
218  /**
219  * @tparam Random the type of the random generator.
220  * @param[in,out] r a random generator.
221  * @param[in] k the index of a digit of the fraction
222  * @return digit number \e k, generating it if necessary.
223  **********************************************************************/
224  template<class Random> unsigned Digit(Random& r, unsigned k) {
225  ExpandTo(r, k + 1);
226  return _f[k];
227  }
228  /**
229  * Add one digit to the fraction.
230  *
231  * @tparam Random the type of the random generator.
232  * @param[in,out] r a random generator.
233  **********************************************************************/
234  template<class Random> void AddDigit(Random& r)
235  { _f.push_back(RandomDigit(r)); }
236  /**
237  * @param[in] k the index of a digit of the fraction
238  * @return a const reference to digit number \e k, without generating new
239  * digits.
240  * @exception std::out_of_range if the digit hasn't been generated.
241  **********************************************************************/
242  const unsigned& RawDigit(unsigned k) const throw()
243  { return (const unsigned&)(_f.at(k)); }
244  /**
245  * @param[in] k the index of a digit of the fraction
246  * @return a non-const reference to digit number \e k, without generating
247  * new digits.
248  * @exception std::out_of_range if the digit hasn't been generated.
249  **********************************************************************/
250  unsigned& RawDigit(unsigned k) throw()
251  { return (unsigned&)(_f.at(k)); }
252  /**
253  * Return to initial state, uniformly distributed in \e n + (0,1). This is
254  * similar to Init but also returns the memory used by the object to the
255  * system. Normally Init should be used.
256  **********************************************************************/
257  void Clear() {
258  std::vector<unsigned> z(0);
259  _n = 0;
260  _s = 1;
261  _f.swap(z);
262  }
263  /**
264  * @return the number of digits in fraction
265  **********************************************************************/
266  unsigned Size() const throw() { return unsigned(_f.size()); }
267  /**
268  * Return the fraction part of the RandomNumber as a floating point number
269  * of type RealType rounded to the nearest multiple of
270  * 1/2<sup><i>p</i></sup>, where \e p =
271  * std::numeric_limits<RealType>::digits, and, if necessary, creating
272  * additional digits of the number.
273  *
274  * @tparam RealType the floating point type to convert to.
275  * @tparam Random the type of the random generator.
276  * @param[in,out] r a random generator for generating the necessary digits.
277  * @return the fraction of the RandomNumber rounded to a RealType.
278  **********************************************************************/
279  template<typename RealType, typename Random> RealType Fraction(Random& r) {
280  STATIC_ASSERT(!std::numeric_limits<RealType>::is_integer,
281  "RandomNumber::Fraction: invalid real type RealType");
282  const int d = std::numeric_limits<RealType>::digits;
283  const int k = (d + bits - 1)/bits;
284  const int kg = (d + bits)/bits; // For guard bit
285  RealType y = 0;
286  if (Digit(r, kg - 1) & (1U << (kg * bits - d - 1)))
287  // if guard bit is set, round up.
288  y += std::pow(RealType(2), -d);
289  const RealType fact = std::pow(RealType(2), -bits);
290  RealType mult = RealType(1);
291  for (int i = 0; i < k; ++i) {
292  mult *= fact;
293  y += mult * RealType(i < k - 1 ? RawDigit(i) :
294  RawDigit(i) & (~0U << (k * bits - d)));
295  }
296  return y;
297  }
298  /**
299  * Return the value of the RandomNumber rounded to nearest floating point
300  * number of type RealType and, if necessary, creating additional digits of
301  * the number.
302  *
303  * @tparam RealType the floating point type to convert to.
304  * @tparam Random the type of the random generator.
305  * @param[in,out] r a random generator for generating the necessary digits.
306  * @return the value of the RandomNumber rounded to a RealType.
307  **********************************************************************/
308  template<typename RealType, class Random> RealType Value(Random& r) {
309  // Ignore the possibility of overflow here (OK because int doesn't
310  // currently overflow any real type). Assume the real type supports
311  // denormalized numbers. Need to treat rounding explicitly since the
312  // missing digits always imply rounding up.
313  STATIC_ASSERT(!std::numeric_limits<RealType>::is_integer,
314  "RandomNumber::Value: invalid real type RealType");
315  const int digits = std::numeric_limits<RealType>::digits,
316  min_exp = std::numeric_limits<RealType>::min_exponent;
317  RealType y;
318  int lead; // Position of leading bit (0.5 = position 0)
319  if (_n) lead = highest_bit_idx(_n);
320  else {
321  int i = 0;
322  while ( Digit(r, i) == 0 && i < (-min_exp)/bits ) ++i;
323  lead = highest_bit_idx(RawDigit(i)) - (i + 1) * bits;
324  // To handle denormalized numbers set lead = max(lead, min_exp)
325  lead = lead > min_exp ? lead : min_exp;
326  }
327  int trail = lead - digits; // Position of guard bit (0.5 = position 0)
328  if (trail > 0) {
329  y = RealType(_n & (~0U << trail));
330  if (_n & (1U << (trail - 1)))
331  y += std::pow(RealType(2), trail);
332  } else {
333  y = RealType(_n);
334  int k = (-trail)/bits; // Byte with guard bit
335  if (Digit(r, k) & (1U << ((k + 1) * bits + trail - 1)))
336  // If guard bit is set, round bit (some subsequent bit will be 1).
337  y += std::pow(RealType(2), trail);
338  // Byte with trailing bit (can be negative)
339  k = (-trail - 1 + bits)/bits - 1;
340  const RealType fact = std::pow(RealType(2), -bits);
341  RealType mult = RealType(1);
342  for (int i = 0; i <= k; ++i) {
343  mult *= fact;
344  y += mult *
345  RealType(i < k ? RawDigit(i) :
346  RawDigit(i) & (~0U << ((k + 1) * bits + trail)));
347  }
348  }
349  if (_s < 0) y *= -1;
350  return y;
351  }
352  /**
353  * Return the range of possible values for the RandomNumber as pair of
354  * doubles. This doesn't create any additional digits of the result and
355  * doesn't try to control roundoff.
356  *
357  * @return a pair denoting the range with first being the lower limit and
358  * second being the upper limit.
359  **********************************************************************/
360  std::pair<double, double> Range() const throw() {
361  double y = _n;
362  const double fact = std::pow(double(2), -bits);
363  double mult = double(1);
364  for (unsigned i = 0; i < Size(); ++i) {
365  mult *= fact;
366  y += mult * RawDigit(i);
367  }
368  return std::pair<double, double>(_s > 0 ? y : -(y + mult),
369  _s > 0 ? (y + mult) : -y);
370  }
371  /**
372  * @tparam Random the type of the random generator.
373  * @param[in,out] r a random generator.
374  * @return a random digit in [0, 2<sup><i>bits</i></sup>).
375  **********************************************************************/
376  template<class Random> static unsigned RandomDigit(Random& r) throw()
377  { return unsigned(r.template Integer<bits>()); }
378 
379  private:
380  /**
381  * The integer part
382  **********************************************************************/
383  unsigned _n;
384  /**
385  * The sign
386  **********************************************************************/
387  int _s;
388  /**
389  * The fraction part
390  **********************************************************************/
391  std::vector<unsigned> _f;
392  /**
393  * Fill RandomNumber to \e k digits.
394  **********************************************************************/
395  template<class Random> void ExpandTo(Random& r, size_t k) {
396  size_t l = _f.size();
397  if (k <= l)
398  return;
399  _f.resize(k);
400  for (size_t i = l; i < k; ++i)
401  _f[i] = RandomDigit(r);
402  }
403  /**
404  * Return index [0..32] of highest bit set. Return 0 if x = 0, 32 is if x
405  * = ~0. (From Algorithms for programmers by Joerg Arndt.)
406  **********************************************************************/
407  static int highest_bit_idx(unsigned x) throw() {
408  if (x == 0) return 0;
409  int r = 1;
410  if (x & 0xffff0000U) { x >>= 16; r += 16; }
411  if (x & 0x0000ff00U) { x >>= 8; r += 8; }
412  if (x & 0x000000f0U) { x >>= 4; r += 4; }
413  if (x & 0x0000000cU) { x >>= 2; r += 2; }
414  if (x & 0x00000002U) { r += 1; }
415  return r;
416  }
417  /**
418  * The number of bits in unsigned.
419  **********************************************************************/
420  static const int w = std::numeric_limits<unsigned>::digits;
421  public:
422  /**
423  * A mask for the digits.
424  **********************************************************************/
425  static const unsigned mask =
426  bits == w ? ~0U : ~(~0U << (bits < w ? bits : 0));
427  };
428 
429  /**
430  * \relates RandomNumber
431  * Print a RandomNumber. Format is n.dddd... where the base for printing is
432  * 2<sup>max(4,<i>b</i>)</sup>. The ... represents an infinite sequence of
433  * ungenerated random digits (uniformly distributed). Thus with \e b = 1,
434  * 0.0... = (0,1/2), 0.00... = (0,1/4), 0.11... = (3/4,1), etc.
435  **********************************************************************/
436  template<int bits>
437  std::ostream& operator<<(std::ostream& os, const RandomNumber<bits>& n) {
438  const std::ios::fmtflags oldflags = os.flags();
439  RandomNumber<bits> t = n;
440  os << (t.Sign() > 0 ? "+" : "-");
441  unsigned i = t.UInteger();
442  os << std::hex << std::setfill('0');
443  if (i == 0)
444  os << "0";
445  else {
446  bool first = true;
447  const int w = std::numeric_limits<unsigned>::digits;
448  const unsigned mask = RandomNumber<bits>::mask;
449  for (int s = ((w + bits - 1)/bits) * bits - bits; s >= 0; s -= bits) {
450  unsigned d = mask & (i >> s);
451  if (d || !first) {
452  if (first) {
453  os << d;
454  first = false;
455  }
456  else
457  os << std::setw((bits+3)/4) << d;
458  }
459  }
460  }
461  os << ".";
462  unsigned s = t.Size();
463  for (unsigned i = 0; i < s; ++i)
464  os << std::setw((bits+3)/4) << t.RawDigit(i);
465  os << "..." << std::setfill(' ');
466  os.flags(oldflags);
467  return os;
468  }
469 
470 } // namespace RandomLib
471 
472 #endif // RANDOMLIB_RANDOMNUMBER_HPP
Header for UniformInteger.
bool LessThan(Random &r, IntType p, IntType q)
Generate random integers, reals, and booleans.
bool GreaterThanEqual(Random &r, IntType p, IntType q)
Infinite precision random numbers.
unsigned UInteger() const
unsigned & RawDigit(unsigned k)
std::pair< double, double > Range() const
const unsigned & RawDigit(unsigned k) const
The partial uniform integer distribution.
void swap(RandomNumber &t)
static unsigned RandomDigit(Random &r)
bool LessThan(Random &r, IntType p0, IntType c, IntType q, UniformInteger< IntType, bits > &j)
#define STATIC_ASSERT(cond, reason)
Definition: MPFRRandom.hpp:31
RealType Value(Random &r)
void swap(RandomLib::RandomCanonical< Generator > &r, RandomLib::RandomCanonical< Generator > &s)
static const unsigned mask
void AddDigit(Random &r)
unsigned Size() const
bool LessThan(Random &r, RandomNumber &t)
unsigned Digit(Random &r, unsigned k)
RealType Fraction(Random &r)
bool GreaterPair(Random &r, RandomNumber &u, RandomNumber &v)