1 // Copyright 2019 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/base/internal/exponential_biased.h"
16
17 #include <stdint.h>
18
19 #include <algorithm>
20 #include <atomic>
21 #include <cmath>
22 #include <limits>
23
24 #include "absl/base/attributes.h"
25 #include "absl/base/optimization.h"
26
27 namespace absl {
28 ABSL_NAMESPACE_BEGIN
29 namespace base_internal {
30
31 // The algorithm generates a random number between 0 and 1 and applies the
32 // inverse cumulative distribution function for an exponential. Specifically:
33 // Let m be the inverse of the sample period, then the probability
34 // distribution function is m*exp(-mx) so the CDF is
35 // p = 1 - exp(-mx), so
36 // q = 1 - p = exp(-mx)
37 // log_e(q) = -mx
38 // -log_e(q)/m = x
39 // log_2(q) * (-log_e(2) * 1/m) = x
40 // In the code, q is actually in the range 1 to 2**26, hence the -26 below
GetSkipCount(int64_t mean)41 int64_t ExponentialBiased::GetSkipCount(int64_t mean) {
42 if (ABSL_PREDICT_FALSE(!initialized_)) {
43 Initialize();
44 }
45
46 uint64_t rng = NextRandom(rng_);
47 rng_ = rng;
48
49 // Take the top 26 bits as the random number
50 // (This plus the 1<<58 sampling bound give a max possible step of
51 // 5194297183973780480 bytes.)
52 // The uint32_t cast is to prevent a (hard-to-reproduce) NAN
53 // under piii debug for some binaries.
54 double q = static_cast<uint32_t>(rng >> (kPrngNumBits - 26)) + 1.0;
55 // Put the computed p-value through the CDF of a geometric.
56 double interval = bias_ + (std::log2(q) - 26) * (-std::log(2.0) * mean);
57 // Very large values of interval overflow int64_t. To avoid that, we will
58 // cheat and clamp any huge values to (int64_t max)/2. This is a potential
59 // source of bias, but the mean would need to be such a large value that it's
60 // not likely to come up. For example, with a mean of 1e18, the probability of
61 // hitting this condition is about 1/1000. For a mean of 1e17, standard
62 // calculators claim that this event won't happen.
63 if (interval > static_cast<double>(std::numeric_limits<int64_t>::max() / 2)) {
64 // Assume huge values are bias neutral, retain bias for next call.
65 return std::numeric_limits<int64_t>::max() / 2;
66 }
67 double value = std::rint(interval);
68 bias_ = interval - value;
69 return value;
70 }
71
GetStride(int64_t mean)72 int64_t ExponentialBiased::GetStride(int64_t mean) {
73 return GetSkipCount(mean - 1) + 1;
74 }
75
Initialize()76 void ExponentialBiased::Initialize() {
77 // We don't get well distributed numbers from `this` so we call NextRandom() a
78 // bunch to mush the bits around. We use a global_rand to handle the case
79 // where the same thread (by memory address) gets created and destroyed
80 // repeatedly.
81 ABSL_CONST_INIT static std::atomic<uint32_t> global_rand(0);
82 uint64_t r = reinterpret_cast<uint64_t>(this) +
83 global_rand.fetch_add(1, std::memory_order_relaxed);
84 for (int i = 0; i < 20; ++i) {
85 r = NextRandom(r);
86 }
87 rng_ = r;
88 initialized_ = true;
89 }
90
91 } // namespace base_internal
92 ABSL_NAMESPACE_END
93 } // namespace absl
94