summaryrefslogtreecommitdiff
path: root/absl/random/gaussian_distribution.h
blob: c1427b06d5496b1fedb6309d9dc4cbb3e6492199 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
// Copyright 2017 The Abseil Authors.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//      https://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#ifndef ABSL_RANDOM_GAUSSIAN_DISTRIBUTION_H_
#define ABSL_RANDOM_GAUSSIAN_DISTRIBUTION_H_

// absl::gaussian_distribution implements the Ziggurat algorithm
// for generating random gaussian numbers.
//
// Implementation based on "The Ziggurat Method for Generating Random Variables"
// by George Marsaglia and Wai Wan Tsang: http://www.jstatsoft.org/v05/i08/
//

#include <cmath>
#include <cstdint>
#include <istream>
#include <limits>
#include <type_traits>

#include "absl/random/internal/fast_uniform_bits.h"
#include "absl/random/internal/generate_real.h"
#include "absl/random/internal/iostream_state_saver.h"

namespace absl {
ABSL_NAMESPACE_BEGIN
namespace random_internal {

// absl::gaussian_distribution_base implements the underlying ziggurat algorithm
// using the ziggurat tables generated by the gaussian_distribution_gentables
// binary.
//
// The specific algorithm has some of the improvements suggested by the
// 2005 paper, "An Improved Ziggurat Method to Generate Normal Random Samples",
// Jurgen A Doornik.  (https://www.doornik.com/research/ziggurat.pdf)
class gaussian_distribution_base {
 public:
  template <typename URBG>
  inline double zignor(URBG& g);  // NOLINT(runtime/references)

 private:
  friend class TableGenerator;

  template <typename URBG>
  inline double zignor_fallback(URBG& g,  // NOLINT(runtime/references)
                                bool neg);

  // Constants used for the gaussian distribution.
  static constexpr double kR = 3.442619855899;  // Start of the tail.
  static constexpr double kRInv = 0.29047645161474317;  // ~= (1.0 / kR) .
  static constexpr double kV = 9.91256303526217e-3;
  static constexpr uint64_t kMask = 0x07f;

  // The ziggurat tables store the pdf(f) and inverse-pdf(x) for equal-area
  // points on one-half of the normal distribution, where the pdf function,
  // pdf = e ^ (-1/2 *x^2), assumes that the mean = 0 & stddev = 1.
  //
  // These tables are just over 2kb in size; larger tables might improve the
  // distributions, but also lead to more cache pollution.
  //
  // x = {3.71308, 3.44261, 3.22308, ..., 0}
  // f = {0.00101, 0.00266, 0.00554, ..., 1}
  struct Tables {
    double x[kMask + 2];
    double f[kMask + 2];
  };
  static const Tables zg_;
  random_internal::FastUniformBits<uint64_t> fast_u64_;
};

}  // namespace random_internal

// absl::gaussian_distribution:
// Generates a number conforming to a Gaussian distribution.
template <typename RealType = double>
class gaussian_distribution : random_internal::gaussian_distribution_base {
 public:
  using result_type = RealType;

  class param_type {
   public:
    using distribution_type = gaussian_distribution;

    explicit param_type(result_type mean = 0, result_type stddev = 1)
        : mean_(mean), stddev_(stddev) {}

    // Returns the mean distribution parameter.  The mean specifies the location
    // of the peak.  The default value is 0.0.
    result_type mean() const { return mean_; }

    // Returns the deviation distribution parameter.  The default value is 1.0.
    result_type stddev() const { return stddev_; }

    friend bool operator==(const param_type& a, const param_type& b) {
      return a.mean_ == b.mean_ && a.stddev_ == b.stddev_;
    }

    friend bool operator!=(const param_type& a, const param_type& b) {
      return !(a == b);
    }

   private:
    result_type mean_;
    result_type stddev_;

    static_assert(
        std::is_floating_point<RealType>::value,
        "Class-template absl::gaussian_distribution<> must be parameterized "
        "using a floating-point type.");
  };

  gaussian_distribution() : gaussian_distribution(0) {}

  explicit gaussian_distribution(result_type mean, result_type stddev = 1)
      : param_(mean, stddev) {}

  explicit gaussian_distribution(const param_type& p) : param_(p) {}

  void reset() {}

  // Generating functions
  template <typename URBG>
  result_type operator()(URBG& g) {  // NOLINT(runtime/references)
    return (*this)(g, param_);
  }

  template <typename URBG>
  result_type operator()(URBG& g,  // NOLINT(runtime/references)
                         const param_type& p);

  param_type param() const { return param_; }
  void param(const param_type& p) { param_ = p; }

  result_type(min)() const {
    return -std::numeric_limits<result_type>::infinity();
  }
  result_type(max)() const {
    return std::numeric_limits<result_type>::infinity();
  }

  result_type mean() const { return param_.mean(); }
  result_type stddev() const { return param_.stddev(); }

  friend bool operator==(const gaussian_distribution& a,
                         const gaussian_distribution& b) {
    return a.param_ == b.param_;
  }
  friend bool operator!=(const gaussian_distribution& a,
                         const gaussian_distribution& b) {
    return a.param_ != b.param_;
  }

 private:
  param_type param_;
};

// --------------------------------------------------------------------------
// Implementation details only below
// --------------------------------------------------------------------------

template <typename RealType>
template <typename URBG>
typename gaussian_distribution<RealType>::result_type
gaussian_distribution<RealType>::operator()(
    URBG& g,  // NOLINT(runtime/references)
    const param_type& p) {
  return p.mean() + p.stddev() * static_cast<result_type>(zignor(g));
}

template <typename CharT, typename Traits, typename RealType>
std::basic_ostream<CharT, Traits>& operator<<(
    std::basic_ostream<CharT, Traits>& os,  // NOLINT(runtime/references)
    const gaussian_distribution<RealType>& x) {
  auto saver = random_internal::make_ostream_state_saver(os);
  os.precision(random_internal::stream_precision_helper<RealType>::kPrecision);
  os << x.mean() << os.fill() << x.stddev();
  return os;
}

template <typename CharT, typename Traits, typename RealType>
std::basic_istream<CharT, Traits>& operator>>(
    std::basic_istream<CharT, Traits>& is,  // NOLINT(runtime/references)
    gaussian_distribution<RealType>& x) {   // NOLINT(runtime/references)
  using result_type = typename gaussian_distribution<RealType>::result_type;
  using param_type = typename gaussian_distribution<RealType>::param_type;

  auto saver = random_internal::make_istream_state_saver(is);
  auto mean = random_internal::read_floating_point<result_type>(is);
  if (is.fail()) return is;
  auto stddev = random_internal::read_floating_point<result_type>(is);
  if (!is.fail()) {
    x.param(param_type(mean, stddev));
  }
  return is;
}

namespace random_internal {

template <typename URBG>
inline double gaussian_distribution_base::zignor_fallback(URBG& g, bool neg) {
  using random_internal::GeneratePositiveTag;
  using random_internal::GenerateRealFromBits;

  // This fallback path happens approximately 0.05% of the time.
  double x, y;
  do {
    // kRInv = 1/r, U(0, 1)
    x = kRInv *
        std::log(GenerateRealFromBits<double, GeneratePositiveTag, false>(
            fast_u64_(g)));
    y = -std::log(
        GenerateRealFromBits<double, GeneratePositiveTag, false>(fast_u64_(g)));
  } while ((y + y) < (x * x));
  return neg ? (x - kR) : (kR - x);
}

template <typename URBG>
inline double gaussian_distribution_base::zignor(
    URBG& g) {  // NOLINT(runtime/references)
  using random_internal::GeneratePositiveTag;
  using random_internal::GenerateRealFromBits;
  using random_internal::GenerateSignedTag;

  while (true) {
    // We use a single uint64_t to generate both a double and a strip.
    // These bits are unused when the generated double is > 1/2^5.
    // This may introduce some bias from the duplicated low bits of small
    // values (those smaller than 1/2^5, which all end up on the left tail).
    uint64_t bits = fast_u64_(g);
    int i = static_cast<int>(bits & kMask);  // pick a random strip
    double j = GenerateRealFromBits<double, GenerateSignedTag, false>(
        bits);  // U(-1, 1)
    const double x = j * zg_.x[i];

    // Retangular box. Handles >97% of all cases.
    // For any given box, this handles between 75% and 99% of values.
    // Equivalent to U(01) < (x[i+1] / x[i]), and when i == 0, ~93.5%
    if (std::abs(x) < zg_.x[i + 1]) {
      return x;
    }

    // i == 0: Base box. Sample using a ratio of uniforms.
    if (i == 0) {
      // This path happens about 0.05% of the time.
      return zignor_fallback(g, j < 0);
    }

    // i > 0: Wedge samples using precomputed values.
    double v = GenerateRealFromBits<double, GeneratePositiveTag, false>(
        fast_u64_(g));  // U(0, 1)
    if ((zg_.f[i + 1] + v * (zg_.f[i] - zg_.f[i + 1])) <
        std::exp(-0.5 * x * x)) {
      return x;
    }

    // The wedge was missed; reject the value and try again.
  }
}

}  // namespace random_internal
ABSL_NAMESPACE_END
}  // namespace absl

#endif  // ABSL_RANDOM_GAUSSIAN_DISTRIBUTION_H_