ExRandom  3.0
unit_normal_dist.hpp
Go to the documentation of this file.
1 /**
2  * @file unit_normal_dist.hpp
3  * @author Charles Karney <charles.karney@sri.com>
4  * @brief Definition of unit_normal_dist
5  *
6  * Copyright (c) Charles Karney (2014-2020) and licensed under the MIT/X11
7  * License. For more information, see http://exrandom.sourceforge.net/
8  */
9 
10 #if !defined(EXRANDOM_UNIT_NORMAL_DIST_HPP)
11 #define EXRANDOM_UNIT_NORMAL_DIST_HPP 1
12 
13 #include <algorithm> // for std::min, std::max
14 
15 #include <exrandom/u_rand.hpp>
16 
17 namespace exrandom {
18 
19  /**
20  * @brief Sample u-rands exactly from the unit normal distribution.
21  *
22  * This samples from the unit normal distribution
23  * P(x) = exp(&minus; x<sup>2</sup>/2) / &radic;(2&pi;).
24  * This implements Algorithm N with improvements due to Du et al. (2020).
25  *
26  * @tparam digit_gen the type of digit generator.
27  *
28  * This class allows a u-rand to be returned via the
29  * unit_normal_dist::generate member function or a floating point result via
30  * the unit_normal_dist::value<RealType> member function.
31  *
32  * See unit_normal_distribution for a simpler interface to sampling normal
33  * deviates. (But, you can't obtain u-rands with this class and the
34  * constructor locks you into a specific floating point type.)
35  */
36  template<typename digit_gen> class unit_normal_dist {
37  public:
38  /**
39  * The constructor.
40  *
41  * @param D a reference to the digit generator to be used.
42  */
43  unit_normal_dist(digit_gen& D) : _D(D), _y(_D), _z(_D), _x(_D) {}
44 
45  /**
46  * Generate the next deviate as a u_rand.
47  *
48  * @tparam Generator the type of g.
49  * @param g the random generator engine.
50  * @param[out] x the u_rand to set.
51  */
52  template<typename Generator>
53  void generate(Generator& g, u_rand<digit_gen>& x) {
54  for (;;) {
55  int k = G(g); // step 1
56  k = S(k); if (k < 0) continue; // step 2
57  x.init(); // step 3
58  int j = k; while (j-- && E(g, x)) {}; // step 4
59  if (!(j < 0)) continue;
60  if (!B(g, x)) continue; // step 5
61  x.set_integer(k); // step 6
62  if (_y.init().less_than_half(g)) x.negate(); // step 7
63  return; // step 8
64  }
65  }
66 
67  /**
68  * Generate the next deviate and round it to a floating point number.
69  *
70  * @tparam RealType the floating point type of the result.
71  * @tparam Generator the type of g.
72  * @param g the random generator engine.
73  * @return the normal deviate.
74  */
75  template<typename RealType, typename Generator>
76  RealType value(Generator& g) {
77  generate(g, _x);
78  return _x.template value<RealType>(g);
79  }
80 
81  /**
82  * @return a reference to the digit generator used in the constructor.
83  */
84  digit_gen& digit_generator() const { return _D; }
85 
86  /**
87  * The base of the digit generator.
88  */
89  static const uint_t base = digit_gen::base;
90  private:
91  static const int bits = digit_gen::bits;
92  static const bool power_of_two = digit_gen::power_of_two;
93  // Disable copy assignment
94  unit_normal_dist& operator=(const unit_normal_dist&);
95  digit_gen& _D;
96  u_rand<digit_gen> _y, _z, _x; // temporary storage
97 
98  // Algorithm H: true with probability exp(-1/2).
99  template<typename Generator>
100  bool H(Generator& g) {
101  if (!_y.init().less_than_half(g)) return true;
102  for (;;) {
103  if (!_z.init().less_than(g, _y)) return false;
104  if (!_y.init().less_than(g, _z)) return true;
105  }
106  }
107 
108  // Step N1: return n >= 0 with prob. exp(-n/2) * (1 - exp(-1/2)).
109  template<typename Generator>
110  int G(Generator& g)
111  { int n = 0; while (H(g)) ++n; return n; }
112 
113  // Step N2: return the square root of n >= 0 if it's a perfect square else -1
114  int S(int n) const {
115  for (int k = 0, k2 = 0; k2 <= n; ++k, k2 += 2*k - 1) {
116  // Here k2 = k * k; note that k^2 - (k - 1)^2 = 2*k - 1
117  if (n == k2) return k;
118  }
119  return -1;
120  }
121 
122  // Algorithm E: true with probability exp(-x) for x in (0, 1).
123  template<typename Generator>
124  bool E(Generator& g, u_rand<digit_gen>& x) {
125  if (!_y.init().less_than(g, x)) return true;
126  for (;;) {
127  if (!_z.init().less_than(g, _y)) return false;
128  if (!_y.init().less_than(g, _z)) return true;
129  }
130  }
131 
132  // Algorithm B: true with prob exp(-x^2 / 2).
133  template<typename Generator>
134  bool B(Generator& g, u_rand<digit_gen>& x) {
135  int n = 0;
136  for (;; ++n) {
137  if (_z.init().less_than_half(g)) break;
138  if (!_z.init().less_than(g, n ? _y : x)) break;
139  if (!_y.init().less_than(g, x)) break;
140  _y.swap(_z); // an efficient way of doing y = z
141  }
142  return (n % 2) == 0;
143  }
144 
145  };
146 
147 }
148 
149 #endif // EXRANDOM_UNIT_NORMAL_DIST_HPP
Sample u-rands exactly from the unit normal distribution.
digit_gen & digit_generator() const
RealType value(Generator &g)
Definition of u_rand.
void generate(Generator &g, u_rand< digit_gen > &x)
void swap(u_rand &t)
Definition: u_rand.hpp:134
u_rand & init()
Definition: u_rand.hpp:123
The common namespace.
Definition: aux_info.hpp:18
The machinery to handle u-rands, arbitrary precision random deviates.
Definition: u_rand.hpp:105
void set_integer(unsigned n)
Definition: u_rand.hpp:155