• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 // Copyright 2017 The Abseil Authors.
2 //
3 // Licensed under the Apache License, Version 2.0 (the "License");
4 // you may not use this file except in compliance with the License.
5 // You may obtain a copy of the License at
6 //
7 //      https://www.apache.org/licenses/LICENSE-2.0
8 //
9 // Unless required by applicable law or agreed to in writing, software
10 // distributed under the License is distributed on an "AS IS" BASIS,
11 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12 // See the License for the specific language governing permissions and
13 // limitations under the License.
14 
15 #ifndef ABSL_RANDOM_INTERNAL_FAST_UNIFORM_BITS_H_
16 #define ABSL_RANDOM_INTERNAL_FAST_UNIFORM_BITS_H_
17 
18 #include <cstddef>
19 #include <cstdint>
20 #include <limits>
21 #include <type_traits>
22 
23 #include "absl/base/config.h"
24 
25 namespace absl {
26 ABSL_NAMESPACE_BEGIN
27 namespace random_internal {
28 // Returns true if the input value is zero or a power of two. Useful for
29 // determining if the range of output values in a URBG
30 template <typename UIntType>
IsPowerOfTwoOrZero(UIntType n)31 constexpr bool IsPowerOfTwoOrZero(UIntType n) {
32   return (n == 0) || ((n & (n - 1)) == 0);
33 }
34 
35 // Computes the length of the range of values producible by the URBG, or returns
36 // zero if that would encompass the entire range of representable values in
37 // URBG::result_type.
38 template <typename URBG>
RangeSize()39 constexpr typename URBG::result_type RangeSize() {
40   using result_type = typename URBG::result_type;
41   return ((URBG::max)() == (std::numeric_limits<result_type>::max)() &&
42           (URBG::min)() == std::numeric_limits<result_type>::lowest())
43              ? result_type{0}
44              : (URBG::max)() - (URBG::min)() + result_type{1};
45 }
46 
47 template <typename UIntType>
LargestPowerOfTwoLessThanOrEqualTo(UIntType n)48 constexpr UIntType LargestPowerOfTwoLessThanOrEqualTo(UIntType n) {
49   return n < 2 ? n : 2 * LargestPowerOfTwoLessThanOrEqualTo(n / 2);
50 }
51 
52 // Given a URBG generating values in the closed interval [Lo, Hi], returns the
53 // largest power of two less than or equal to `Hi - Lo + 1`.
54 template <typename URBG>
PowerOfTwoSubRangeSize()55 constexpr typename URBG::result_type PowerOfTwoSubRangeSize() {
56   return LargestPowerOfTwoLessThanOrEqualTo(RangeSize<URBG>());
57 }
58 
59 // Computes the floor of the log. (i.e., std::floor(std::log2(N));
60 template <typename UIntType>
IntegerLog2(UIntType n)61 constexpr UIntType IntegerLog2(UIntType n) {
62   return (n <= 1) ? 0 : 1 + IntegerLog2(n / 2);
63 }
64 
65 // Returns the number of bits of randomness returned through
66 // `PowerOfTwoVariate(urbg)`.
67 template <typename URBG>
NumBits()68 constexpr size_t NumBits() {
69   return RangeSize<URBG>() == 0
70              ? std::numeric_limits<typename URBG::result_type>::digits
71              : IntegerLog2(PowerOfTwoSubRangeSize<URBG>());
72 }
73 
74 // Given a shift value `n`, constructs a mask with exactly the low `n` bits set.
75 // If `n == 0`, all bits are set.
76 template <typename UIntType>
MaskFromShift(UIntType n)77 constexpr UIntType MaskFromShift(UIntType n) {
78   return ((n % std::numeric_limits<UIntType>::digits) == 0)
79              ? ~UIntType{0}
80              : (UIntType{1} << n) - UIntType{1};
81 }
82 
83 // FastUniformBits implements a fast path to acquire uniform independent bits
84 // from a type which conforms to the [rand.req.urbg] concept.
85 // Parameterized by:
86 //  `UIntType`: the result (output) type
87 //
88 // The std::independent_bits_engine [rand.adapt.ibits] adaptor can be
89 // instantiated from an existing generator through a copy or a move. It does
90 // not, however, facilitate the production of pseudorandom bits from an un-owned
91 // generator that will outlive the std::independent_bits_engine instance.
92 template <typename UIntType = uint64_t>
93 class FastUniformBits {
94  public:
95   using result_type = UIntType;
96 
result_type(min)97   static constexpr result_type(min)() { return 0; }
result_type(max)98   static constexpr result_type(max)() {
99     return (std::numeric_limits<result_type>::max)();
100   }
101 
102   template <typename URBG>
103   result_type operator()(URBG& g);  // NOLINT(runtime/references)
104 
105  private:
106   static_assert(std::is_unsigned<UIntType>::value,
107                 "Class-template FastUniformBits<> must be parameterized using "
108                 "an unsigned type.");
109 
110   // PowerOfTwoVariate() generates a single random variate, always returning a
111   // value in the half-open interval `[0, PowerOfTwoSubRangeSize<URBG>())`. If
112   // the URBG already generates values in a power-of-two range, the generator
113   // itself is used. Otherwise, we use rejection sampling on the largest
114   // possible power-of-two-sized subrange.
115   struct PowerOfTwoTag {};
116   struct RejectionSamplingTag {};
117   template <typename URBG>
PowerOfTwoVariate(URBG & g)118   static typename URBG::result_type PowerOfTwoVariate(
119       URBG& g) {  // NOLINT(runtime/references)
120     using tag =
121         typename std::conditional<IsPowerOfTwoOrZero(RangeSize<URBG>()),
122                                   PowerOfTwoTag, RejectionSamplingTag>::type;
123     return PowerOfTwoVariate(g, tag{});
124   }
125 
126   template <typename URBG>
PowerOfTwoVariate(URBG & g,PowerOfTwoTag)127   static typename URBG::result_type PowerOfTwoVariate(
128       URBG& g,  // NOLINT(runtime/references)
129       PowerOfTwoTag) {
130     return g() - (URBG::min)();
131   }
132 
133   template <typename URBG>
PowerOfTwoVariate(URBG & g,RejectionSamplingTag)134   static typename URBG::result_type PowerOfTwoVariate(
135       URBG& g,  // NOLINT(runtime/references)
136       RejectionSamplingTag) {
137     // Use rejection sampling to ensure uniformity across the range.
138     typename URBG::result_type u;
139     do {
140       u = g() - (URBG::min)();
141     } while (u >= PowerOfTwoSubRangeSize<URBG>());
142     return u;
143   }
144 
145   // Generate() generates a random value, dispatched on whether
146   // the underlying URBG must loop over multiple calls or not.
147   template <typename URBG>
148   result_type Generate(URBG& g,  // NOLINT(runtime/references)
149                        std::true_type /* avoid_looping */);
150 
151   template <typename URBG>
152   result_type Generate(URBG& g,  // NOLINT(runtime/references)
153                        std::false_type /* avoid_looping */);
154 };
155 
156 template <typename UIntType>
157 template <typename URBG>
158 typename FastUniformBits<UIntType>::result_type
operator()159 FastUniformBits<UIntType>::operator()(URBG& g) {  // NOLINT(runtime/references)
160   // kRangeMask is the mask used when sampling variates from the URBG when the
161   // width of the URBG range is not a power of 2.
162   // Y = (2 ^ kRange) - 1
163   static_assert((URBG::max)() > (URBG::min)(),
164                 "URBG::max and URBG::min may not be equal.");
165   using urbg_result_type = typename URBG::result_type;
166   constexpr urbg_result_type kRangeMask =
167       RangeSize<URBG>() == 0
168           ? (std::numeric_limits<urbg_result_type>::max)()
169           : static_cast<urbg_result_type>(PowerOfTwoSubRangeSize<URBG>() - 1);
170   return Generate(g, std::integral_constant<bool, (kRangeMask >= (max)())>{});
171 }
172 
173 template <typename UIntType>
174 template <typename URBG>
175 typename FastUniformBits<UIntType>::result_type
Generate(URBG & g,std::true_type)176 FastUniformBits<UIntType>::Generate(URBG& g,  // NOLINT(runtime/references)
177                                     std::true_type /* avoid_looping */) {
178   // The width of the result_type is less than than the width of the random bits
179   // provided by URBG.  Thus, generate a single value and then simply mask off
180   // the required bits.
181 
182   return PowerOfTwoVariate(g) & (max)();
183 }
184 
185 template <typename UIntType>
186 template <typename URBG>
187 typename FastUniformBits<UIntType>::result_type
Generate(URBG & g,std::false_type)188 FastUniformBits<UIntType>::Generate(URBG& g,  // NOLINT(runtime/references)
189                                     std::false_type /* avoid_looping */) {
190   // See [rand.adapt.ibits] for more details on the constants calculated below.
191   //
192   // It is preferable to use roughly the same number of bits from each generator
193   // call, however this is only possible when the number of bits provided by the
194   // URBG is a divisor of the number of bits in `result_type`. In all other
195   // cases, the number of bits used cannot always be the same, but it can be
196   // guaranteed to be off by at most 1. Thus we run two loops, one with a
197   // smaller bit-width size (`kSmallWidth`) and one with a larger width size
198   // (satisfying `kLargeWidth == kSmallWidth + 1`). The loops are run
199   // `kSmallIters` and `kLargeIters` times respectively such
200   // that
201   //
202   //    `kTotalWidth == kSmallIters * kSmallWidth
203   //                    + kLargeIters * kLargeWidth`
204   //
205   // where `kTotalWidth` is the total number of bits in `result_type`.
206   //
207   constexpr size_t kTotalWidth = std::numeric_limits<result_type>::digits;
208   constexpr size_t kUrbgWidth = NumBits<URBG>();
209   constexpr size_t kTotalIters =
210       kTotalWidth / kUrbgWidth + (kTotalWidth % kUrbgWidth != 0);
211   constexpr size_t kSmallWidth = kTotalWidth / kTotalIters;
212   constexpr size_t kLargeWidth = kSmallWidth + 1;
213   //
214   // Because `kLargeWidth == kSmallWidth + 1`, it follows that
215   //
216   //     `kTotalWidth == kTotalIters * kSmallWidth + kLargeIters`
217   //
218   // and therefore
219   //
220   //     `kLargeIters == kTotalWidth % kSmallWidth`
221   //
222   // Intuitively, each iteration with the large width accounts for one unit
223   // of the remainder when `kTotalWidth` is divided by `kSmallWidth`. As
224   // mentioned above, if the URBG width is a divisor of `kTotalWidth`, then
225   // there would be no need for any large iterations (i.e., one loop would
226   // suffice), and indeed, in this case, `kLargeIters` would be zero.
227   constexpr size_t kLargeIters = kTotalWidth % kSmallWidth;
228   constexpr size_t kSmallIters =
229       (kTotalWidth - (kLargeWidth * kLargeIters)) / kSmallWidth;
230 
231   static_assert(
232       kTotalWidth == kSmallIters * kSmallWidth + kLargeIters * kLargeWidth,
233       "Error in looping constant calculations.");
234 
235   result_type s = 0;
236 
237   constexpr size_t kSmallShift = kSmallWidth % kTotalWidth;
238   constexpr result_type kSmallMask = MaskFromShift(result_type{kSmallShift});
239   for (size_t n = 0; n < kSmallIters; ++n) {
240     s = (s << kSmallShift) +
241         (static_cast<result_type>(PowerOfTwoVariate(g)) & kSmallMask);
242   }
243 
244   constexpr size_t kLargeShift = kLargeWidth % kTotalWidth;
245   constexpr result_type kLargeMask = MaskFromShift(result_type{kLargeShift});
246   for (size_t n = 0; n < kLargeIters; ++n) {
247     s = (s << kLargeShift) +
248         (static_cast<result_type>(PowerOfTwoVariate(g)) & kLargeMask);
249   }
250 
251   static_assert(
252       kLargeShift == kSmallShift + 1 ||
253           (kLargeShift == 0 &&
254            kSmallShift == std::numeric_limits<result_type>::digits - 1),
255       "Error in looping constant calculations");
256 
257   return s;
258 }
259 
260 }  // namespace random_internal
261 ABSL_NAMESPACE_END
262 }  // namespace absl
263 
264 #endif  // ABSL_RANDOM_INTERNAL_FAST_UNIFORM_BITS_H_
265