diff options
Diffstat (limited to 'third_party/abseil_cpp/absl/random/gaussian_distribution.h')
| -rw-r--r-- | third_party/abseil_cpp/absl/random/gaussian_distribution.h | 275 |
1 files changed, 0 insertions, 275 deletions
diff --git a/third_party/abseil_cpp/absl/random/gaussian_distribution.h b/third_party/abseil_cpp/absl/random/gaussian_distribution.h deleted file mode 100644 index 4b07a5c0af9d..000000000000 --- a/third_party/abseil_cpp/absl/random/gaussian_distribution.h +++ /dev/null @@ -1,275 +0,0 @@ -// 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/base/config.h" -#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 ABSL_DLL 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_ |