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 #include "absl/random/uniform_real_distribution.h"
16
17 #include <cfloat>
18 #include <cmath>
19 #include <cstdint>
20 #include <iterator>
21 #include <random>
22 #include <sstream>
23 #include <string>
24 #include <type_traits>
25 #include <vector>
26
27 #include "gmock/gmock.h"
28 #include "gtest/gtest.h"
29 #include "absl/log/log.h"
30 #include "absl/numeric/internal/representation.h"
31 #include "absl/random/internal/chi_square.h"
32 #include "absl/random/internal/distribution_test_util.h"
33 #include "absl/random/internal/pcg_engine.h"
34 #include "absl/random/internal/sequence_urbg.h"
35 #include "absl/random/random.h"
36 #include "absl/strings/str_cat.h"
37
38 // NOTES:
39 // * Some documentation on generating random real values suggests that
40 // it is possible to use std::nextafter(b, DBL_MAX) to generate a value on
41 // the closed range [a, b]. Unfortunately, that technique is not universally
42 // reliable due to floating point quantization.
43 //
44 // * absl::uniform_real_distribution<float> generates between 2^28 and 2^29
45 // distinct floating point values in the range [0, 1).
46 //
47 // * absl::uniform_real_distribution<float> generates at least 2^23 distinct
48 // floating point values in the range [1, 2). This should be the same as
49 // any other range covered by a single exponent in IEEE 754.
50 //
51 // * absl::uniform_real_distribution<double> generates more than 2^52 distinct
52 // values in the range [0, 1), and should generate at least 2^52 distinct
53 // values in the range of [1, 2).
54 //
55
56 namespace {
57
58 template <typename RealType>
59 class UniformRealDistributionTest : public ::testing::Test {};
60
61 // double-double arithmetic is not supported well by either GCC or Clang; see
62 // https://gcc.gnu.org/bugzilla/show_bug.cgi?id=99048,
63 // https://bugs.llvm.org/show_bug.cgi?id=49131, and
64 // https://bugs.llvm.org/show_bug.cgi?id=49132. Don't bother running these tests
65 // with double doubles until compiler support is better.
66 using RealTypes =
67 std::conditional<absl::numeric_internal::IsDoubleDouble(),
68 ::testing::Types<float, double>,
69 ::testing::Types<float, double, long double>>::type;
70
71 TYPED_TEST_SUITE(UniformRealDistributionTest, RealTypes);
72
TYPED_TEST(UniformRealDistributionTest,ParamSerializeTest)73 TYPED_TEST(UniformRealDistributionTest, ParamSerializeTest) {
74 #if (defined(__i386__) || defined(_M_IX86)) && FLT_EVAL_METHOD != 0
75 // We're using an x87-compatible FPU, and intermediate operations are
76 // performed with 80-bit floats. This produces slightly different results from
77 // what we expect below.
78 GTEST_SKIP()
79 << "Skipping the test because we detected x87 floating-point semantics";
80 #endif
81 using DistributionType = absl::uniform_real_distribution<TypeParam>;
82 using real_type = TypeParam;
83 using param_type = typename DistributionType::param_type;
84
85 constexpr const real_type kMax = std::numeric_limits<real_type>::max();
86 constexpr const real_type kMin = std::numeric_limits<real_type>::min();
87 constexpr const real_type kEpsilon =
88 std::numeric_limits<real_type>::epsilon();
89 constexpr const real_type kLowest =
90 std::numeric_limits<real_type>::lowest(); // -max
91
92 const real_type kDenormMax = std::nextafter(kMin, real_type{0});
93 const real_type kOneMinusE =
94 std::nextafter(real_type{1}, real_type{0}); // 1 - epsilon
95
96 constexpr const real_type kTwo60{1152921504606846976}; // 2^60
97
98 constexpr int kCount = 1000;
99 absl::InsecureBitGen gen;
100 for (const auto& param : {
101 param_type(),
102 param_type(real_type{0}, real_type{1}),
103 param_type(real_type(-0.1), real_type(0.1)),
104 param_type(real_type(0.05), real_type(0.12)),
105 param_type(real_type(-0.05), real_type(0.13)),
106 param_type(real_type(-0.05), real_type(-0.02)),
107 // range = 0
108 param_type(real_type(2.0), real_type(2.0)), // Same
109 // double range = 0
110 // 2^60 , 2^60 + 2^6
111 param_type(kTwo60, real_type(1152921504606847040)),
112 // 2^60 , 2^60 + 2^7
113 param_type(kTwo60, real_type(1152921504606847104)),
114 // double range = 2^8
115 // 2^60 , 2^60 + 2^8
116 param_type(kTwo60, real_type(1152921504606847232)),
117 // float range = 0
118 // 2^60 , 2^60 + 2^36
119 param_type(kTwo60, real_type(1152921573326323712)),
120 // 2^60 , 2^60 + 2^37
121 param_type(kTwo60, real_type(1152921642045800448)),
122 // float range = 2^38
123 // 2^60 , 2^60 + 2^38
124 param_type(kTwo60, real_type(1152921779484753920)),
125 // Limits
126 param_type(0, kMax),
127 param_type(kLowest, 0),
128 param_type(0, kMin),
129 param_type(0, kEpsilon),
130 param_type(-kEpsilon, kEpsilon),
131 param_type(0, kOneMinusE),
132 param_type(0, kDenormMax),
133 }) {
134 // Validate parameters.
135 const auto a = param.a();
136 const auto b = param.b();
137 DistributionType before(a, b);
138 EXPECT_EQ(before.a(), param.a());
139 EXPECT_EQ(before.b(), param.b());
140
141 {
142 DistributionType via_param(param);
143 EXPECT_EQ(via_param, before);
144 }
145
146 std::stringstream ss;
147 ss << before;
148 DistributionType after(real_type(1.0), real_type(3.1));
149
150 EXPECT_NE(before.a(), after.a());
151 EXPECT_NE(before.b(), after.b());
152 EXPECT_NE(before.param(), after.param());
153 EXPECT_NE(before, after);
154
155 ss >> after;
156
157 EXPECT_EQ(before.a(), after.a());
158 EXPECT_EQ(before.b(), after.b());
159 EXPECT_EQ(before.param(), after.param());
160 EXPECT_EQ(before, after);
161
162 // Smoke test.
163 auto sample_min = after.max();
164 auto sample_max = after.min();
165 for (int i = 0; i < kCount; i++) {
166 auto sample = after(gen);
167 // Failure here indicates a bug in uniform_real_distribution::operator(),
168 // or bad parameters--range too large, etc.
169 if (after.min() == after.max()) {
170 EXPECT_EQ(sample, after.min());
171 } else {
172 EXPECT_GE(sample, after.min());
173 EXPECT_LT(sample, after.max());
174 }
175 if (sample > sample_max) {
176 sample_max = sample;
177 }
178 if (sample < sample_min) {
179 sample_min = sample;
180 }
181 }
182
183 if (!std::is_same<real_type, long double>::value) {
184 // static_cast<double>(long double) can overflow.
185 LOG(INFO) << "Range: " << static_cast<double>(sample_min) << ", "
186 << static_cast<double>(sample_max);
187 }
188 }
189 }
190
191 #ifdef _MSC_VER
192 #pragma warning(push)
193 #pragma warning(disable:4756) // Constant arithmetic overflow.
194 #endif
TYPED_TEST(UniformRealDistributionTest,ViolatesPreconditionsDeathTest)195 TYPED_TEST(UniformRealDistributionTest, ViolatesPreconditionsDeathTest) {
196 using DistributionType = absl::uniform_real_distribution<TypeParam>;
197 using real_type = TypeParam;
198
199 #if GTEST_HAS_DEATH_TEST
200 // Hi < Lo
201 EXPECT_DEBUG_DEATH({ DistributionType dist(10.0, 1.0); }, "");
202
203 // Hi - Lo > numeric_limits<>::max()
204 EXPECT_DEBUG_DEATH(
205 {
206 DistributionType dist(std::numeric_limits<real_type>::lowest(),
207 std::numeric_limits<real_type>::max());
208 },
209 "");
210
211 // kEpsilon guarantees that max + kEpsilon = inf.
212 const auto kEpsilon = std::nexttoward(
213 (std::numeric_limits<real_type>::max() -
214 std::nexttoward(std::numeric_limits<real_type>::max(), 0.0)) /
215 2,
216 std::numeric_limits<real_type>::max());
217 EXPECT_DEBUG_DEATH(
218 {
219 DistributionType dist(-kEpsilon, std::numeric_limits<real_type>::max());
220 },
221 "");
222 EXPECT_DEBUG_DEATH(
223 {
224 DistributionType dist(std::numeric_limits<real_type>::lowest(),
225 kEpsilon);
226 },
227 "");
228
229 #endif // GTEST_HAS_DEATH_TEST
230 #if defined(NDEBUG)
231 // opt-mode, for invalid parameters, will generate a garbage value,
232 // but should not enter an infinite loop.
233 absl::InsecureBitGen gen;
234 {
235 DistributionType dist(10.0, 1.0);
236 auto x = dist(gen);
237 EXPECT_FALSE(std::isnan(x)) << x;
238 }
239 {
240 DistributionType dist(std::numeric_limits<real_type>::lowest(),
241 std::numeric_limits<real_type>::max());
242 auto x = dist(gen);
243 // Infinite result.
244 EXPECT_FALSE(std::isfinite(x)) << x;
245 }
246 #endif // NDEBUG
247 }
248 #ifdef _MSC_VER
249 #pragma warning(pop) // warning(disable:4756)
250 #endif
251
TYPED_TEST(UniformRealDistributionTest,TestMoments)252 TYPED_TEST(UniformRealDistributionTest, TestMoments) {
253 using DistributionType = absl::uniform_real_distribution<TypeParam>;
254
255 constexpr int kSize = 1000000;
256 std::vector<double> values(kSize);
257
258 // We use a fixed bit generator for distribution accuracy tests. This allows
259 // these tests to be deterministic, while still testing the qualify of the
260 // implementation.
261 absl::random_internal::pcg64_2018_engine rng{0x2B7E151628AED2A6};
262
263 DistributionType dist;
264 for (int i = 0; i < kSize; i++) {
265 values[i] = dist(rng);
266 }
267
268 const auto moments =
269 absl::random_internal::ComputeDistributionMoments(values);
270 EXPECT_NEAR(0.5, moments.mean, 0.01);
271 EXPECT_NEAR(1 / 12.0, moments.variance, 0.015);
272 EXPECT_NEAR(0.0, moments.skewness, 0.02);
273 EXPECT_NEAR(9 / 5.0, moments.kurtosis, 0.015);
274 }
275
TYPED_TEST(UniformRealDistributionTest,ChiSquaredTest50)276 TYPED_TEST(UniformRealDistributionTest, ChiSquaredTest50) {
277 using DistributionType = absl::uniform_real_distribution<TypeParam>;
278 using param_type = typename DistributionType::param_type;
279
280 using absl::random_internal::kChiSquared;
281
282 constexpr size_t kTrials = 100000;
283 constexpr int kBuckets = 50;
284 constexpr double kExpected =
285 static_cast<double>(kTrials) / static_cast<double>(kBuckets);
286
287 // 1-in-100000 threshold, but remember, there are about 8 tests
288 // in this file. And the test could fail for other reasons.
289 // Empirically validated with --runs_per_test=10000.
290 const int kThreshold =
291 absl::random_internal::ChiSquareValue(kBuckets - 1, 0.999999);
292
293 // We use a fixed bit generator for distribution accuracy tests. This allows
294 // these tests to be deterministic, while still testing the qualify of the
295 // implementation.
296 absl::random_internal::pcg64_2018_engine rng{0x2B7E151628AED2A6};
297
298 for (const auto& param : {param_type(0, 1), param_type(5, 12),
299 param_type(-5, 13), param_type(-5, -2)}) {
300 const double min_val = param.a();
301 const double max_val = param.b();
302 const double factor = kBuckets / (max_val - min_val);
303
304 std::vector<int32_t> counts(kBuckets, 0);
305 DistributionType dist(param);
306 for (size_t i = 0; i < kTrials; i++) {
307 auto x = dist(rng);
308 auto bucket = static_cast<size_t>((x - min_val) * factor);
309 counts[bucket]++;
310 }
311
312 double chi_square = absl::random_internal::ChiSquareWithExpected(
313 std::begin(counts), std::end(counts), kExpected);
314 if (chi_square > kThreshold) {
315 double p_value =
316 absl::random_internal::ChiSquarePValue(chi_square, kBuckets);
317
318 // Chi-squared test failed. Output does not appear to be uniform.
319 std::string msg;
320 for (const auto& a : counts) {
321 absl::StrAppend(&msg, a, "\n");
322 }
323 absl::StrAppend(&msg, kChiSquared, " p-value ", p_value, "\n");
324 absl::StrAppend(&msg, "High ", kChiSquared, " value: ", chi_square, " > ",
325 kThreshold);
326 LOG(INFO) << msg;
327 FAIL() << msg;
328 }
329 }
330 }
331
TYPED_TEST(UniformRealDistributionTest,StabilityTest)332 TYPED_TEST(UniformRealDistributionTest, StabilityTest) {
333 using DistributionType = absl::uniform_real_distribution<TypeParam>;
334 using real_type = TypeParam;
335
336 // absl::uniform_real_distribution stability relies only on
337 // random_internal::GenerateRealFromBits.
338 absl::random_internal::sequence_urbg urbg(
339 {0x0003eb76f6f7f755ull, 0xFFCEA50FDB2F953Bull, 0xC332DDEFBE6C5AA5ull,
340 0x6558218568AB9702ull, 0x2AEF7DAD5B6E2F84ull, 0x1521B62829076170ull,
341 0xECDD4775619F1510ull, 0x13CCA830EB61BD96ull, 0x0334FE1EAA0363CFull,
342 0xB5735C904C70A239ull, 0xD59E9E0BCBAADE14ull, 0xEECC86BC60622CA7ull});
343
344 std::vector<int> output(12);
345
346 DistributionType dist;
347 std::generate(std::begin(output), std::end(output), [&] {
348 return static_cast<int>(real_type(1000000) * dist(urbg));
349 });
350
351 EXPECT_THAT(
352 output, //
353 testing::ElementsAre(59, 999246, 762494, 395876, 167716, 82545, 925251,
354 77341, 12527, 708791, 834451, 932808));
355 }
356
TEST(UniformRealDistributionTest,AlgorithmBounds)357 TEST(UniformRealDistributionTest, AlgorithmBounds) {
358 absl::uniform_real_distribution<double> dist;
359
360 {
361 // This returns the smallest value >0 from absl::uniform_real_distribution.
362 absl::random_internal::sequence_urbg urbg({0x0000000000000001ull});
363 double a = dist(urbg);
364 EXPECT_EQ(a, 5.42101086242752217004e-20);
365 }
366
367 {
368 // This returns a value very near 0.5 from absl::uniform_real_distribution.
369 absl::random_internal::sequence_urbg urbg({0x7fffffffffffffefull});
370 double a = dist(urbg);
371 EXPECT_EQ(a, 0.499999999999999944489);
372 }
373 {
374 // This returns a value very near 0.5 from absl::uniform_real_distribution.
375 absl::random_internal::sequence_urbg urbg({0x8000000000000000ull});
376 double a = dist(urbg);
377 EXPECT_EQ(a, 0.5);
378 }
379
380 {
381 // This returns the largest value <1 from absl::uniform_real_distribution.
382 absl::random_internal::sequence_urbg urbg({0xFFFFFFFFFFFFFFEFull});
383 double a = dist(urbg);
384 EXPECT_EQ(a, 0.999999999999999888978);
385 }
386 {
387 // This *ALSO* returns the largest value <1.
388 absl::random_internal::sequence_urbg urbg({0xFFFFFFFFFFFFFFFFull});
389 double a = dist(urbg);
390 EXPECT_EQ(a, 0.999999999999999888978);
391 }
392 }
393
394 } // namespace
395