1 | // Class template uniform_int_distribution -*- C++ -*- |
2 | |
3 | // Copyright (C) 2009-2021 Free Software Foundation, Inc. |
4 | // |
5 | // This file is part of the GNU ISO C++ Library. This library is free |
6 | // software; you can redistribute it and/or modify it under the |
7 | // terms of the GNU General Public License as published by the |
8 | // Free Software Foundation; either version 3, or (at your option) |
9 | // any later version. |
10 | |
11 | // This library is distributed in the hope that it will be useful, |
12 | // but WITHOUT ANY WARRANTY; without even the implied warranty of |
13 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
14 | // GNU General Public License for more details. |
15 | |
16 | // Under Section 7 of GPL version 3, you are granted additional |
17 | // permissions described in the GCC Runtime Library Exception, version |
18 | // 3.1, as published by the Free Software Foundation. |
19 | |
20 | // You should have received a copy of the GNU General Public License and |
21 | // a copy of the GCC Runtime Library Exception along with this program; |
22 | // see the files COPYING3 and COPYING.RUNTIME respectively. If not, see |
23 | // <http://www.gnu.org/licenses/>. |
24 | |
25 | /** |
26 | * @file bits/uniform_int_dist.h |
27 | * This is an internal header file, included by other library headers. |
28 | * Do not attempt to use it directly. @headername{random} |
29 | */ |
30 | |
31 | #ifndef _GLIBCXX_BITS_UNIFORM_INT_DIST_H |
32 | #define _GLIBCXX_BITS_UNIFORM_INT_DIST_H |
33 | |
34 | #include <type_traits> |
35 | #include <ext/numeric_traits.h> |
36 | #if __cplusplus > 201703L |
37 | # include <concepts> |
38 | #endif |
39 | #include <bits/concept_check.h> // __glibcxx_function_requires |
40 | |
41 | namespace std _GLIBCXX_VISIBILITY(default) |
42 | { |
43 | _GLIBCXX_BEGIN_NAMESPACE_VERSION |
44 | |
45 | #ifdef __cpp_lib_concepts |
46 | /// Requirements for a uniform random bit generator. |
47 | template<typename _Gen> |
48 | concept uniform_random_bit_generator |
49 | = invocable<_Gen&> && unsigned_integral<invoke_result_t<_Gen&>> |
50 | && requires |
51 | { |
52 | { _Gen::min() } -> same_as<invoke_result_t<_Gen&>>; |
53 | { _Gen::max() } -> same_as<invoke_result_t<_Gen&>>; |
54 | requires bool_constant<(_Gen::min() < _Gen::max())>::value; |
55 | }; |
56 | #endif |
57 | |
58 | namespace __detail |
59 | { |
60 | // Determine whether number is a power of two. |
61 | // This is true for zero, which is OK because we want _Power_of_2(n+1) |
62 | // to be true if n==numeric_limits<_Tp>::max() and so n+1 wraps around. |
63 | template<typename _Tp> |
64 | constexpr bool |
65 | _Power_of_2(_Tp __x) |
66 | { |
67 | return ((__x - 1) & __x) == 0; |
68 | } |
69 | } |
70 | |
71 | /** |
72 | * @brief Uniform discrete distribution for random numbers. |
73 | * A discrete random distribution on the range @f$[min, max]@f$ with equal |
74 | * probability throughout the range. |
75 | */ |
76 | template<typename _IntType = int> |
77 | class uniform_int_distribution |
78 | { |
79 | static_assert(std::is_integral<_IntType>::value, |
80 | "template argument must be an integral type" ); |
81 | |
82 | public: |
83 | /** The type of the range of the distribution. */ |
84 | typedef _IntType result_type; |
85 | /** Parameter type. */ |
86 | struct param_type |
87 | { |
88 | typedef uniform_int_distribution<_IntType> distribution_type; |
89 | |
90 | param_type() : param_type(0) { } |
91 | |
92 | explicit |
93 | param_type(_IntType __a, |
94 | _IntType __b = __gnu_cxx::__int_traits<_IntType>::__max) |
95 | : _M_a(__a), _M_b(__b) |
96 | { |
97 | __glibcxx_assert(_M_a <= _M_b); |
98 | } |
99 | |
100 | result_type |
101 | a() const |
102 | { return _M_a; } |
103 | |
104 | result_type |
105 | b() const |
106 | { return _M_b; } |
107 | |
108 | friend bool |
109 | operator==(const param_type& __p1, const param_type& __p2) |
110 | { return __p1._M_a == __p2._M_a && __p1._M_b == __p2._M_b; } |
111 | |
112 | friend bool |
113 | operator!=(const param_type& __p1, const param_type& __p2) |
114 | { return !(__p1 == __p2); } |
115 | |
116 | private: |
117 | _IntType _M_a; |
118 | _IntType _M_b; |
119 | }; |
120 | |
121 | public: |
122 | /** |
123 | * @brief Constructs a uniform distribution object. |
124 | */ |
125 | uniform_int_distribution() : uniform_int_distribution(0) { } |
126 | |
127 | /** |
128 | * @brief Constructs a uniform distribution object. |
129 | */ |
130 | explicit |
131 | uniform_int_distribution(_IntType __a, |
132 | _IntType __b |
133 | = __gnu_cxx::__int_traits<_IntType>::__max) |
134 | : _M_param(__a, __b) |
135 | { } |
136 | |
137 | explicit |
138 | uniform_int_distribution(const param_type& __p) |
139 | : _M_param(__p) |
140 | { } |
141 | |
142 | /** |
143 | * @brief Resets the distribution state. |
144 | * |
145 | * Does nothing for the uniform integer distribution. |
146 | */ |
147 | void |
148 | reset() { } |
149 | |
150 | result_type |
151 | a() const |
152 | { return _M_param.a(); } |
153 | |
154 | result_type |
155 | b() const |
156 | { return _M_param.b(); } |
157 | |
158 | /** |
159 | * @brief Returns the parameter set of the distribution. |
160 | */ |
161 | param_type |
162 | param() const |
163 | { return _M_param; } |
164 | |
165 | /** |
166 | * @brief Sets the parameter set of the distribution. |
167 | * @param __param The new parameter set of the distribution. |
168 | */ |
169 | void |
170 | param(const param_type& __param) |
171 | { _M_param = __param; } |
172 | |
173 | /** |
174 | * @brief Returns the inclusive lower bound of the distribution range. |
175 | */ |
176 | result_type |
177 | min() const |
178 | { return this->a(); } |
179 | |
180 | /** |
181 | * @brief Returns the inclusive upper bound of the distribution range. |
182 | */ |
183 | result_type |
184 | max() const |
185 | { return this->b(); } |
186 | |
187 | /** |
188 | * @brief Generating functions. |
189 | */ |
190 | template<typename _UniformRandomBitGenerator> |
191 | result_type |
192 | operator()(_UniformRandomBitGenerator& __urng) |
193 | { return this->operator()(__urng, _M_param); } |
194 | |
195 | template<typename _UniformRandomBitGenerator> |
196 | result_type |
197 | operator()(_UniformRandomBitGenerator& __urng, |
198 | const param_type& __p); |
199 | |
200 | template<typename _ForwardIterator, |
201 | typename _UniformRandomBitGenerator> |
202 | void |
203 | __generate(_ForwardIterator __f, _ForwardIterator __t, |
204 | _UniformRandomBitGenerator& __urng) |
205 | { this->__generate(__f, __t, __urng, _M_param); } |
206 | |
207 | template<typename _ForwardIterator, |
208 | typename _UniformRandomBitGenerator> |
209 | void |
210 | __generate(_ForwardIterator __f, _ForwardIterator __t, |
211 | _UniformRandomBitGenerator& __urng, |
212 | const param_type& __p) |
213 | { this->__generate_impl(__f, __t, __urng, __p); } |
214 | |
215 | template<typename _UniformRandomBitGenerator> |
216 | void |
217 | __generate(result_type* __f, result_type* __t, |
218 | _UniformRandomBitGenerator& __urng, |
219 | const param_type& __p) |
220 | { this->__generate_impl(__f, __t, __urng, __p); } |
221 | |
222 | /** |
223 | * @brief Return true if two uniform integer distributions have |
224 | * the same parameters. |
225 | */ |
226 | friend bool |
227 | operator==(const uniform_int_distribution& __d1, |
228 | const uniform_int_distribution& __d2) |
229 | { return __d1._M_param == __d2._M_param; } |
230 | |
231 | private: |
232 | template<typename _ForwardIterator, |
233 | typename _UniformRandomBitGenerator> |
234 | void |
235 | __generate_impl(_ForwardIterator __f, _ForwardIterator __t, |
236 | _UniformRandomBitGenerator& __urng, |
237 | const param_type& __p); |
238 | |
239 | param_type _M_param; |
240 | |
241 | // Lemire's nearly divisionless algorithm. |
242 | // Returns an unbiased random number from __g downscaled to [0,__range) |
243 | // using an unsigned type _Wp twice as wide as unsigned type _Up. |
244 | template<typename _Wp, typename _Urbg, typename _Up> |
245 | static _Up |
246 | _S_nd(_Urbg& __g, _Up __range) |
247 | { |
248 | using _Up_traits = __gnu_cxx::__int_traits<_Up>; |
249 | using _Wp_traits = __gnu_cxx::__int_traits<_Wp>; |
250 | static_assert(!_Up_traits::__is_signed, "U must be unsigned" ); |
251 | static_assert(!_Wp_traits::__is_signed, "W must be unsigned" ); |
252 | static_assert(_Wp_traits::__digits == (2 * _Up_traits::__digits), |
253 | "W must be twice as wide as U" ); |
254 | |
255 | // reference: Fast Random Integer Generation in an Interval |
256 | // ACM Transactions on Modeling and Computer Simulation 29 (1), 2019 |
257 | // https://arxiv.org/abs/1805.10941 |
258 | _Wp __product = _Wp(__g()) * _Wp(__range); |
259 | _Up __low = _Up(__product); |
260 | if (__low < __range) |
261 | { |
262 | _Up __threshold = -__range % __range; |
263 | while (__low < __threshold) |
264 | { |
265 | __product = _Wp(__g()) * _Wp(__range); |
266 | __low = _Up(__product); |
267 | } |
268 | } |
269 | return __product >> _Up_traits::__digits; |
270 | } |
271 | }; |
272 | |
273 | template<typename _IntType> |
274 | template<typename _UniformRandomBitGenerator> |
275 | typename uniform_int_distribution<_IntType>::result_type |
276 | uniform_int_distribution<_IntType>:: |
277 | operator()(_UniformRandomBitGenerator& __urng, |
278 | const param_type& __param) |
279 | { |
280 | typedef typename _UniformRandomBitGenerator::result_type _Gresult_type; |
281 | typedef typename make_unsigned<result_type>::type __utype; |
282 | typedef typename common_type<_Gresult_type, __utype>::type __uctype; |
283 | |
284 | constexpr __uctype __urngmin = _UniformRandomBitGenerator::min(); |
285 | constexpr __uctype __urngmax = _UniformRandomBitGenerator::max(); |
286 | static_assert( __urngmin < __urngmax, |
287 | "Uniform random bit generator must define min() < max()" ); |
288 | constexpr __uctype __urngrange = __urngmax - __urngmin; |
289 | |
290 | const __uctype __urange |
291 | = __uctype(__param.b()) - __uctype(__param.a()); |
292 | |
293 | __uctype __ret; |
294 | if (__urngrange > __urange) |
295 | { |
296 | // downscaling |
297 | |
298 | const __uctype __uerange = __urange + 1; // __urange can be zero |
299 | |
300 | #if defined __UINT64_TYPE__ && defined __UINT32_TYPE__ |
301 | #if __SIZEOF_INT128__ |
302 | if _GLIBCXX17_CONSTEXPR (__urngrange == __UINT64_MAX__) |
303 | { |
304 | // __urng produces values that use exactly 64-bits, |
305 | // so use 128-bit integers to downscale to desired range. |
306 | __UINT64_TYPE__ __u64erange = __uerange; |
307 | __ret = _S_nd<unsigned __int128>(__urng, __u64erange); |
308 | } |
309 | else |
310 | #endif |
311 | if _GLIBCXX17_CONSTEXPR (__urngrange == __UINT32_MAX__) |
312 | { |
313 | // __urng produces values that use exactly 32-bits, |
314 | // so use 64-bit integers to downscale to desired range. |
315 | __UINT32_TYPE__ __u32erange = __uerange; |
316 | __ret = _S_nd<__UINT64_TYPE__>(__urng, __u32erange); |
317 | } |
318 | else |
319 | #endif |
320 | { |
321 | // fallback case (2 divisions) |
322 | const __uctype __scaling = __urngrange / __uerange; |
323 | const __uctype __past = __uerange * __scaling; |
324 | do |
325 | __ret = __uctype(__urng()) - __urngmin; |
326 | while (__ret >= __past); |
327 | __ret /= __scaling; |
328 | } |
329 | } |
330 | else if (__urngrange < __urange) |
331 | { |
332 | // upscaling |
333 | /* |
334 | Note that every value in [0, urange] |
335 | can be written uniquely as |
336 | |
337 | (urngrange + 1) * high + low |
338 | |
339 | where |
340 | |
341 | high in [0, urange / (urngrange + 1)] |
342 | |
343 | and |
344 | |
345 | low in [0, urngrange]. |
346 | */ |
347 | __uctype __tmp; // wraparound control |
348 | do |
349 | { |
350 | const __uctype __uerngrange = __urngrange + 1; |
351 | __tmp = (__uerngrange * operator() |
352 | (__urng, param_type(0, __urange / __uerngrange))); |
353 | __ret = __tmp + (__uctype(__urng()) - __urngmin); |
354 | } |
355 | while (__ret > __urange || __ret < __tmp); |
356 | } |
357 | else |
358 | __ret = __uctype(__urng()) - __urngmin; |
359 | |
360 | return __ret + __param.a(); |
361 | } |
362 | |
363 | |
364 | template<typename _IntType> |
365 | template<typename _ForwardIterator, |
366 | typename _UniformRandomBitGenerator> |
367 | void |
368 | uniform_int_distribution<_IntType>:: |
369 | __generate_impl(_ForwardIterator __f, _ForwardIterator __t, |
370 | _UniformRandomBitGenerator& __urng, |
371 | const param_type& __param) |
372 | { |
373 | __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) |
374 | typedef typename _UniformRandomBitGenerator::result_type _Gresult_type; |
375 | typedef typename make_unsigned<result_type>::type __utype; |
376 | typedef typename common_type<_Gresult_type, __utype>::type __uctype; |
377 | |
378 | static_assert( __urng.min() < __urng.max(), |
379 | "Uniform random bit generator must define min() < max()" ); |
380 | |
381 | constexpr __uctype __urngmin = __urng.min(); |
382 | constexpr __uctype __urngmax = __urng.max(); |
383 | constexpr __uctype __urngrange = __urngmax - __urngmin; |
384 | const __uctype __urange |
385 | = __uctype(__param.b()) - __uctype(__param.a()); |
386 | |
387 | __uctype __ret; |
388 | |
389 | if (__urngrange > __urange) |
390 | { |
391 | if (__detail::_Power_of_2(__urngrange + 1) |
392 | && __detail::_Power_of_2(__urange + 1)) |
393 | { |
394 | while (__f != __t) |
395 | { |
396 | __ret = __uctype(__urng()) - __urngmin; |
397 | *__f++ = (__ret & __urange) + __param.a(); |
398 | } |
399 | } |
400 | else |
401 | { |
402 | // downscaling |
403 | const __uctype __uerange = __urange + 1; // __urange can be zero |
404 | const __uctype __scaling = __urngrange / __uerange; |
405 | const __uctype __past = __uerange * __scaling; |
406 | while (__f != __t) |
407 | { |
408 | do |
409 | __ret = __uctype(__urng()) - __urngmin; |
410 | while (__ret >= __past); |
411 | *__f++ = __ret / __scaling + __param.a(); |
412 | } |
413 | } |
414 | } |
415 | else if (__urngrange < __urange) |
416 | { |
417 | // upscaling |
418 | /* |
419 | Note that every value in [0, urange] |
420 | can be written uniquely as |
421 | |
422 | (urngrange + 1) * high + low |
423 | |
424 | where |
425 | |
426 | high in [0, urange / (urngrange + 1)] |
427 | |
428 | and |
429 | |
430 | low in [0, urngrange]. |
431 | */ |
432 | __uctype __tmp; // wraparound control |
433 | while (__f != __t) |
434 | { |
435 | do |
436 | { |
437 | constexpr __uctype __uerngrange = __urngrange + 1; |
438 | __tmp = (__uerngrange * operator() |
439 | (__urng, param_type(0, __urange / __uerngrange))); |
440 | __ret = __tmp + (__uctype(__urng()) - __urngmin); |
441 | } |
442 | while (__ret > __urange || __ret < __tmp); |
443 | *__f++ = __ret; |
444 | } |
445 | } |
446 | else |
447 | while (__f != __t) |
448 | *__f++ = __uctype(__urng()) - __urngmin + __param.a(); |
449 | } |
450 | |
451 | // operator!= and operator<< and operator>> are defined in <bits/random.h> |
452 | |
453 | _GLIBCXX_END_NAMESPACE_VERSION |
454 | } // namespace std |
455 | |
456 | #endif |
457 | |