blob: 4eb6f754eb8080c6fe652ef0b0fffb3d9fa0e5cd [file] [log] [blame]
terelius84e78f92015-12-10 09:50:551/*
2 * Copyright (c) 2015 The WebRTC project authors. All Rights Reserved.
3 *
4 * Use of this source code is governed by a BSD-style license
5 * that can be found in the LICENSE file in the root of the source
6 * tree. An additional intellectual property rights grant can be found
7 * in the file PATENTS. All contributing project authors may
8 * be found in the AUTHORS file in the root of the source tree.
9 */
10
Jonas Olssona4d87372019-07-05 17:08:3311#include "rtc_base/random.h"
12
terelius84e78f92015-12-10 09:50:5513#include <math.h>
14
15#include <limits>
16#include <vector>
17
Steve Anton10542f22019-01-11 17:11:0018#include "rtc_base/numerics/math_utils.h" // unsigned difference
Mirko Bonadei92ea95e2017-09-15 04:47:3119#include "test/gtest.h"
terelius84e78f92015-12-10 09:50:5520
21namespace webrtc {
22
23namespace {
24// Computes the positive remainder of x/n.
25template <typename T>
26T fdiv_remainder(T x, T n) {
kwiberg352444f2016-11-28 23:58:5327 RTC_CHECK_GE(n, 0);
terelius84e78f92015-12-10 09:50:5528 T remainder = x % n;
29 if (remainder < 0)
30 remainder += n;
31 return remainder;
32}
33} // namespace
34
35// Sample a number of random integers of type T. Divide them into buckets
36// based on the remainder when dividing by bucket_count and check that each
37// bucket gets roughly the expected number of elements.
38template <typename T>
39void UniformBucketTest(T bucket_count, int samples, Random* prng) {
40 std::vector<int> buckets(bucket_count, 0);
41
42 uint64_t total_values = 1ull << (std::numeric_limits<T>::digits +
43 std::numeric_limits<T>::is_signed);
44 T upper_limit =
45 std::numeric_limits<T>::max() -
46 static_cast<T>(total_values % static_cast<uint64_t>(bucket_count));
47 ASSERT_GT(upper_limit, std::numeric_limits<T>::max() / 2);
48
49 for (int i = 0; i < samples; i++) {
50 T sample;
51 do {
52 // We exclude a few numbers from the range so that it is divisible by
53 // the number of buckets. If we are unlucky and hit one of the excluded
54 // numbers we just resample. Note that if the number of buckets is a
55 // power of 2, then we don't have to exclude anything.
56 sample = prng->Rand<T>();
57 } while (sample > upper_limit);
58 buckets[fdiv_remainder(sample, bucket_count)]++;
59 }
60
61 for (T i = 0; i < bucket_count; i++) {
62 // Expect the result to be within 3 standard deviations of the mean.
63 EXPECT_NEAR(buckets[i], samples / bucket_count,
64 3 * sqrt(samples / bucket_count));
65 }
66}
67
68TEST(RandomNumberGeneratorTest, BucketTestSignedChar) {
69 Random prng(7297352569824ull);
70 UniformBucketTest<signed char>(64, 640000, &prng);
71 UniformBucketTest<signed char>(11, 440000, &prng);
72 UniformBucketTest<signed char>(3, 270000, &prng);
73}
74
75TEST(RandomNumberGeneratorTest, BucketTestUnsignedChar) {
76 Random prng(7297352569824ull);
77 UniformBucketTest<unsigned char>(64, 640000, &prng);
78 UniformBucketTest<unsigned char>(11, 440000, &prng);
79 UniformBucketTest<unsigned char>(3, 270000, &prng);
80}
81
82TEST(RandomNumberGeneratorTest, BucketTestSignedShort) {
83 Random prng(7297352569824ull);
84 UniformBucketTest<int16_t>(64, 640000, &prng);
85 UniformBucketTest<int16_t>(11, 440000, &prng);
86 UniformBucketTest<int16_t>(3, 270000, &prng);
87}
88
89TEST(RandomNumberGeneratorTest, BucketTestUnsignedShort) {
90 Random prng(7297352569824ull);
91 UniformBucketTest<uint16_t>(64, 640000, &prng);
92 UniformBucketTest<uint16_t>(11, 440000, &prng);
93 UniformBucketTest<uint16_t>(3, 270000, &prng);
94}
95
96TEST(RandomNumberGeneratorTest, BucketTestSignedInt) {
97 Random prng(7297352569824ull);
98 UniformBucketTest<signed int>(64, 640000, &prng);
99 UniformBucketTest<signed int>(11, 440000, &prng);
100 UniformBucketTest<signed int>(3, 270000, &prng);
101}
102
103TEST(RandomNumberGeneratorTest, BucketTestUnsignedInt) {
104 Random prng(7297352569824ull);
105 UniformBucketTest<unsigned int>(64, 640000, &prng);
106 UniformBucketTest<unsigned int>(11, 440000, &prng);
107 UniformBucketTest<unsigned int>(3, 270000, &prng);
108}
109
110// The range of the random numbers is divided into bucket_count intervals
111// of consecutive numbers. Check that approximately equally many numbers
112// from each inteval are generated.
113void BucketTestSignedInterval(unsigned int bucket_count,
114 unsigned int samples,
115 int32_t low,
116 int32_t high,
117 int sigma_level,
118 Random* prng) {
119 std::vector<unsigned int> buckets(bucket_count, 0);
120
121 ASSERT_GE(high, low);
122 ASSERT_GE(bucket_count, 2u);
Artem Titov9d777622020-09-18 16:23:08123 uint32_t interval = webrtc_impl::unsigned_difference<int32_t>(high, low) + 1;
terelius84e78f92015-12-10 09:50:55124 uint32_t numbers_per_bucket;
125 if (interval == 0) {
126 // The computation high - low + 1 should be 2^32 but overflowed
127 // Hence, bucket_count must be a power of 2
128 ASSERT_EQ(bucket_count & (bucket_count - 1), 0u);
129 numbers_per_bucket = (0x80000000u / bucket_count) * 2;
130 } else {
131 ASSERT_EQ(interval % bucket_count, 0u);
132 numbers_per_bucket = interval / bucket_count;
133 }
134
135 for (unsigned int i = 0; i < samples; i++) {
136 int32_t sample = prng->Rand(low, high);
137 EXPECT_LE(low, sample);
138 EXPECT_GE(high, sample);
Artem Titov9d777622020-09-18 16:23:08139 buckets[webrtc_impl::unsigned_difference<int32_t>(sample, low) /
140 numbers_per_bucket]++;
terelius84e78f92015-12-10 09:50:55141 }
142
143 for (unsigned int i = 0; i < bucket_count; i++) {
144 // Expect the result to be within 3 standard deviations of the mean,
145 // or more generally, within sigma_level standard deviations of the mean.
146 double mean = static_cast<double>(samples) / bucket_count;
147 EXPECT_NEAR(buckets[i], mean, sigma_level * sqrt(mean));
148 }
149}
150
151// The range of the random numbers is divided into bucket_count intervals
152// of consecutive numbers. Check that approximately equally many numbers
153// from each inteval are generated.
154void BucketTestUnsignedInterval(unsigned int bucket_count,
155 unsigned int samples,
156 uint32_t low,
157 uint32_t high,
158 int sigma_level,
159 Random* prng) {
160 std::vector<unsigned int> buckets(bucket_count, 0);
161
162 ASSERT_GE(high, low);
163 ASSERT_GE(bucket_count, 2u);
tereliusd802b5b2016-03-01 19:07:34164 uint32_t interval = high - low + 1;
terelius84e78f92015-12-10 09:50:55165 uint32_t numbers_per_bucket;
166 if (interval == 0) {
167 // The computation high - low + 1 should be 2^32 but overflowed
168 // Hence, bucket_count must be a power of 2
169 ASSERT_EQ(bucket_count & (bucket_count - 1), 0u);
170 numbers_per_bucket = (0x80000000u / bucket_count) * 2;
171 } else {
172 ASSERT_EQ(interval % bucket_count, 0u);
173 numbers_per_bucket = interval / bucket_count;
174 }
175
176 for (unsigned int i = 0; i < samples; i++) {
177 uint32_t sample = prng->Rand(low, high);
178 EXPECT_LE(low, sample);
179 EXPECT_GE(high, sample);
tereliusd802b5b2016-03-01 19:07:34180 buckets[(sample - low) / numbers_per_bucket]++;
terelius84e78f92015-12-10 09:50:55181 }
182
183 for (unsigned int i = 0; i < bucket_count; i++) {
184 // Expect the result to be within 3 standard deviations of the mean,
185 // or more generally, within sigma_level standard deviations of the mean.
186 double mean = static_cast<double>(samples) / bucket_count;
187 EXPECT_NEAR(buckets[i], mean, sigma_level * sqrt(mean));
188 }
189}
190
191TEST(RandomNumberGeneratorTest, UniformUnsignedInterval) {
192 Random prng(299792458ull);
193 BucketTestUnsignedInterval(2, 100000, 0, 1, 3, &prng);
194 BucketTestUnsignedInterval(7, 100000, 1, 14, 3, &prng);
195 BucketTestUnsignedInterval(11, 100000, 1000, 1010, 3, &prng);
196 BucketTestUnsignedInterval(100, 100000, 0, 99, 3, &prng);
197 BucketTestUnsignedInterval(2, 100000, 0, 4294967295, 3, &prng);
198 BucketTestUnsignedInterval(17, 100000, 455, 2147484110, 3, &prng);
199 // 99.7% of all samples will be within 3 standard deviations of the mean,
200 // but since we test 1000 buckets we allow an interval of 4 sigma.
201 BucketTestUnsignedInterval(1000, 1000000, 0, 2147483999, 4, &prng);
202}
203
oprypin8ad0e582017-09-05 10:00:37204TEST(RandomNumberGeneratorTest, UniformSignedInterval) {
terelius84e78f92015-12-10 09:50:55205 Random prng(66260695729ull);
206 BucketTestSignedInterval(2, 100000, 0, 1, 3, &prng);
207 BucketTestSignedInterval(7, 100000, -2, 4, 3, &prng);
208 BucketTestSignedInterval(11, 100000, 1000, 1010, 3, &prng);
209 BucketTestSignedInterval(100, 100000, 0, 99, 3, &prng);
210 BucketTestSignedInterval(2, 100000, std::numeric_limits<int32_t>::min(),
211 std::numeric_limits<int32_t>::max(), 3, &prng);
212 BucketTestSignedInterval(17, 100000, -1073741826, 1073741829, 3, &prng);
213 // 99.7% of all samples will be within 3 standard deviations of the mean,
214 // but since we test 1000 buckets we allow an interval of 4 sigma.
215 BucketTestSignedInterval(1000, 1000000, -352, 2147483647, 4, &prng);
216}
217
218// The range of the random numbers is divided into bucket_count intervals
219// of consecutive numbers. Check that approximately equally many numbers
220// from each inteval are generated.
221void BucketTestFloat(unsigned int bucket_count,
222 unsigned int samples,
223 int sigma_level,
224 Random* prng) {
225 ASSERT_GE(bucket_count, 2u);
226 std::vector<unsigned int> buckets(bucket_count, 0);
227
228 for (unsigned int i = 0; i < samples; i++) {
229 uint32_t sample = bucket_count * prng->Rand<float>();
230 EXPECT_LE(0u, sample);
231 EXPECT_GE(bucket_count - 1, sample);
232 buckets[sample]++;
233 }
234
235 for (unsigned int i = 0; i < bucket_count; i++) {
236 // Expect the result to be within 3 standard deviations of the mean,
237 // or more generally, within sigma_level standard deviations of the mean.
238 double mean = static_cast<double>(samples) / bucket_count;
239 EXPECT_NEAR(buckets[i], mean, sigma_level * sqrt(mean));
240 }
241}
242
243TEST(RandomNumberGeneratorTest, UniformFloatInterval) {
244 Random prng(1380648813ull);
245 BucketTestFloat(100, 100000, 3, &prng);
246 // 99.7% of all samples will be within 3 standard deviations of the mean,
247 // but since we test 1000 buckets we allow an interval of 4 sigma.
248 // BucketTestSignedInterval(1000, 1000000, -352, 2147483647, 4, &prng);
249}
250
251TEST(RandomNumberGeneratorTest, SignedHasSameBitPattern) {
252 Random prng_signed(66738480ull), prng_unsigned(66738480ull);
253
254 for (int i = 0; i < 1000; i++) {
255 signed int s = prng_signed.Rand<signed int>();
256 unsigned int u = prng_unsigned.Rand<unsigned int>();
257 EXPECT_EQ(u, static_cast<unsigned int>(s));
258 }
259
260 for (int i = 0; i < 1000; i++) {
261 int16_t s = prng_signed.Rand<int16_t>();
262 uint16_t u = prng_unsigned.Rand<uint16_t>();
263 EXPECT_EQ(u, static_cast<uint16_t>(s));
264 }
265
266 for (int i = 0; i < 1000; i++) {
267 signed char s = prng_signed.Rand<signed char>();
268 unsigned char u = prng_unsigned.Rand<unsigned char>();
269 EXPECT_EQ(u, static_cast<unsigned char>(s));
270 }
271}
272
273TEST(RandomNumberGeneratorTest, Gaussian) {
274 const int kN = 100000;
275 const int kBuckets = 100;
276 const double kMean = 49;
277 const double kStddev = 10;
278
279 Random prng(1256637061);
280
281 std::vector<unsigned int> buckets(kBuckets, 0);
282 for (int i = 0; i < kN; i++) {
283 int index = prng.Gaussian(kMean, kStddev) + 0.5;
284 if (index >= 0 && index < kBuckets) {
285 buckets[index]++;
286 }
287 }
288
289 const double kPi = 3.14159265358979323846;
290 const double kScale = 1 / (kStddev * sqrt(2.0 * kPi));
291 const double kDiv = -2.0 * kStddev * kStddev;
292 for (int n = 0; n < kBuckets; ++n) {
293 // Use Simpsons rule to estimate the probability that a random gaussian
294 // sample is in the interval [n-0.5, n+0.5].
295 double f_left = kScale * exp((n - kMean - 0.5) * (n - kMean - 0.5) / kDiv);
296 double f_mid = kScale * exp((n - kMean) * (n - kMean) / kDiv);
297 double f_right = kScale * exp((n - kMean + 0.5) * (n - kMean + 0.5) / kDiv);
298 double normal_dist = (f_left + 4 * f_mid + f_right) / 6;
299 // Expect the number of samples to be within 3 standard deviations
300 // (rounded up) of the expected number of samples in the bucket.
301 EXPECT_NEAR(buckets[n], kN * normal_dist, 3 * sqrt(kN * normal_dist) + 1);
302 }
303}
304
305} // namespace webrtc