aarand
Aaron's random distribution functions
Loading...
Searching...
No Matches
aarand.hpp
Go to the documentation of this file.
1#ifndef AARAND_AARAND_HPP
2#define AARAND_AARAND_HPP
3
4#include <cmath>
5#include <limits>
6#include <stdexcept>
7#include <algorithm>
8#include <numeric>
9#include <type_traits>
10
20namespace aarand {
21
33template<typename Output_ = double, class Engine_>
34Output_ standard_uniform(Engine_& eng) {
35 typedef typename Engine_::result_type R;
36 static_assert(std::numeric_limits<R>::is_integer, "RNG engine must yield integer results");
37
38 // Can't be bothered to figure out whether the range fits into 'R' for signed values.
39 // So instead, we just require unsigned integers, where the range will always fit.
40 static_assert(!std::numeric_limits<R>::is_signed, "RNG engine must yield unsigned integers");
41
42 // Stolen from Boost, see https://www.boost.org/doc/libs/1_67_0/boost/random/uniform_01.hpp
43 // The +1 probably doesn't matter for 64-bit generators, but is helpful for engines with
44 // fewer output bits, to reduce the (small) probability of sampling 1's.
45 constexpr Output_ ONE = 1;
46 constexpr Output_ factor = ONE / (static_cast<Output_>(Engine_::max() - Engine_::min()) + ONE);
47
48 // Note that it still might be possible to get a result = 1, depending on
49 // the numerical precision used to compute the product; hence the loop.
50 Output_ result;
51 do {
52 result = static_cast<Output_>(eng() - Engine_::min()) * factor;
53 } while (result == ONE);
54
55 return result;
56}
57
61// Some of the functions below log-transform a uniform random variable.
62// However, standard_uniform() has a small chance of returning zero, resulting in an undesirable -Inf after log.
63// To avoid this, any time we sample zero, we roll again.
64// Note to self: don't try to do something like replacing zeros with 1 to get the log-transform to work,
65// as this introduces inflated occurrences of log(1); better to just reroll to get a true sampling from (0, 1).
66template<typename Output_, class Engine_>
67Output_ non_zero_uniform(Engine_& eng) {
68 Output_ val;
69 do {
71 } while (val == 0);
72 return val;
73}
88template<typename Output_ = double, class Engine_>
89std::pair<Output_, Output_> standard_normal(Engine_& eng) {
90 constexpr Output_ PI = 3.14159265358979323846;
91 constexpr Output_ TWO = 2;
92
93 // Box-Muller gives us two random values at a time.
94 const Output_ constant = std::sqrt(-TWO * std::log(non_zero_uniform<Output_>(eng)));
95 const Output_ angle = TWO * PI * standard_uniform<Output_>(eng);
96 return std::make_pair(constant * std::sin(angle), constant * std::cos(angle));
97}
98
110template<typename Output_ = double, class Engine_>
111Output_ standard_exponential(Engine_& eng) {
112 return -std::log(non_zero_uniform<Output_>(eng));
113}
114
125template<typename Output_ = int, class Engine_>
126Output_ discrete_uniform(Engine_& eng, const Output_ bound) {
127 static_assert(std::numeric_limits<Output_>::is_integer);
128 if (bound <= 0) {
129 throw std::runtime_error("'bound' should be a positive integer");
130 }
131
132 typedef typename Engine_::result_type R;
133 static_assert(std::numeric_limits<R>::is_integer);
134 static_assert(!std::numeric_limits<R>::is_signed); // don't want to figure out how to store the range if it might not fit into R.
135
136 constexpr R range = Engine_::max() - Engine_::min();
137 if (static_cast<typename std::make_unsigned<Output_>::type>(bound) > range) { // force an unsigned comparison.
138 throw std::runtime_error("'bound' should be less than the RNG range");
139 }
140
141 R draw = eng() - Engine_::min();
142
143 // Conservative shortcut to avoid an extra modulo operation in computing
144 // 'limit' if 'draw' is below 'limit'. This is based on the observation
145 // that 'range - bound <= limit', so any condition that triggers the loop
146 // will also pass this check. Allows early return when 'range >> bound'.
147 if (draw > range - bound) {
148
149 // The limit is necessary to provide uniformity in the presence of the
150 // modulus. The idea is to re-sample if we get a draw above the limit.
151 // Technically this can have problems as bound approaches range, in which
152 // case we might end up discarding a lot of the sample space... but this
153 // is unlikely to happen in practice, and even if it does, it's a rejection
154 // rate that's guaranteed to be less than 50%, so whatever.
155 //
156 // Note that the +1 is necessary because range is inclusive but bound is not.
157 const R limit = range - ((range % bound) + 1);
158
159 // In addition, we don't have to deal with the crap about combining draws
160 // to get enough entropy, which is 90% of the Boost implementation.
161 while (draw > limit) {
162 draw = eng() - Engine_::min();
163 }
164 }
165
166 return draw % bound;
167}
168
180template<class InputIterator_, typename Length_, class Engine_>
181void shuffle(const InputIterator_ values, const Length_ n, Engine_& eng) {
182 if (n <= 1) {
183 return;
184 }
185
186 const Length_ last = n - 1;
187 for (Length_ i = 0; i < last; ++i) {
188 const auto chosen = discrete_uniform(eng, n - i);
189 if (chosen) {
190 const auto current = values + i;
191 using std::swap;
192 swap(*current, *(current + chosen));
193 }
194 }
195}
196
212template<class InputIterator_, typename Length_, class OutputIterator_, class Engine_>
213void sample(InputIterator_ values, const Length_ n, const Length_ s, OutputIterator_ output, Engine_& eng) {
214 if (!s) {
215 return;
216 }
217
218 auto remaining = s;
219 for (Length_ i = 0; i < n; ++i, ++values) {
220 const Length_ denom = n - i;
221 const double threshold = static_cast<double>(remaining) / denom;
222 if (threshold >= 1) {
223 // Once remaining >= denom, all remaining values must be selected.
224 // Both values will drop at the same rate so threshold will always be >= 1 in subsequent loops.
225 std::copy_n(values, denom, output);
226 return;
227 }
228
229 if (standard_uniform(eng) <= threshold) {
230 *output = *values;
231 ++output;
232 --remaining;
233 if (!remaining) {
234 return;
235 }
236 }
237 }
238}
239
254template<typename Length_, class OutputIterator_, class Engine_>
255void sample(const Length_ bound, const Length_ s, OutputIterator_ output, Engine_& eng) {
256 if (!s) {
257 return;
258 }
259
260 auto remaining = s;
261 for (Length_ i = 0; i < bound; ++i) {
262 const Length_ denom = bound - i;
263 const double threshold = static_cast<double>(remaining) / denom;
264 if (threshold >= 1) {
265 // Once remaining >= denom, all remaining indices must be selected.
266 // Both values will drop at the same rate so threshold will always be >= 1 in subsequent loops.
267 std::iota(output, output + denom, i);
268 return;
269 }
270
271 if (standard_uniform(eng) <= threshold) {
272 *output = i;
273 ++output;
274 --remaining;
275 if (!remaining) {
276 return;
277 }
278 }
279 }
280}
281
282}
283
284#endif
Aaron's random distribution functions.
void sample(InputIterator_ values, const Length_ n, const Length_ s, OutputIterator_ output, Engine_ &eng)
Definition aarand.hpp:213
void shuffle(const InputIterator_ values, const Length_ n, Engine_ &eng)
Definition aarand.hpp:181
std::pair< Output_, Output_ > standard_normal(Engine_ &eng)
Definition aarand.hpp:89
Output_ standard_exponential(Engine_ &eng)
Definition aarand.hpp:111
Output_ discrete_uniform(Engine_ &eng, const Output_ bound)
Definition aarand.hpp:126
Output_ standard_uniform(Engine_ &eng)
Definition aarand.hpp:34