absl/random/gaussian_distribution.h

0001 // Copyright 2017 The Abseil Authors.
0002 //
0003 // Licensed under the Apache License, Version 2.0 (the "License");
0004 // you may not use this file except in compliance with the License.
0005 // You may obtain a copy of the License at
0006 //
0007 //      https://www.apache.org/licenses/LICENSE-2.0
0008 //
0009 // Unless required by applicable law or agreed to in writing, software
0010 // distributed under the License is distributed on an "AS IS" BASIS,
0011 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
0012 // See the License for the specific language governing permissions and
0013 // limitations under the License.
0014
0015 #ifndef ABSL_RANDOM_GAUSSIAN_DISTRIBUTION_H_
0016 #define ABSL_RANDOM_GAUSSIAN_DISTRIBUTION_H_
0017
0018 // absl::gaussian_distribution implements the Ziggurat algorithm
0019 // for generating random gaussian numbers.
0020 //
0021 // Implementation based on "The Ziggurat Method for Generating Random Variables"
0022 // by George Marsaglia and Wai Wan Tsang: http://www.jstatsoft.org/v05/i08/
0023 //
0024
0025 #include <cmath>
0026 #include <cstdint>
0027 #include <istream>
0028 #include <limits>
0029 #include <type_traits>
0030
0031 #include "absl/base/config.h"
0032 #include "absl/random/internal/fast_uniform_bits.h"
0033 #include "absl/random/internal/generate_real.h"
0034 #include "absl/random/internal/iostream_state_saver.h"
0035
0036 namespace absl {
0037 ABSL_NAMESPACE_BEGIN
0038 namespace random_internal {
0039
0040 // absl::gaussian_distribution_base implements the underlying ziggurat algorithm
0041 // using the ziggurat tables generated by the gaussian_distribution_gentables
0042 // binary.
0043 //
0044 // The specific algorithm has some of the improvements suggested by the
0045 // 2005 paper, "An Improved Ziggurat Method to Generate Normal Random Samples",
0046 // Jurgen A Doornik.  (https://www.doornik.com/research/ziggurat.pdf)
0047 class ABSL_DLL gaussian_distribution_base {
0048  public:
0049   template <typename URBG>
0050   inline double zignor(URBG& g);  // NOLINT(runtime/references)
0051
0052  private:
0053   friend class TableGenerator;
0054
0055   template <typename URBG>
0056   inline double zignor_fallback(URBG& g,  // NOLINT(runtime/references)
0057                                 bool neg);
0058
0059   // Constants used for the gaussian distribution.
0060   static constexpr double kR = 3.442619855899;  // Start of the tail.
0061   static constexpr double kRInv = 0.29047645161474317;  // ~= (1.0 / kR) .
0062   static constexpr double kV = 9.91256303526217e-3;
0063   static constexpr uint64_t kMask = 0x07f;
0064
0065   // The ziggurat tables store the pdf(f) and inverse-pdf(x) for equal-area
0066   // points on one-half of the normal distribution, where the pdf function,
0067   // pdf = e ^ (-1/2 *x^2), assumes that the mean = 0 & stddev = 1.
0068   //
0069   // These tables are just over 2kb in size; larger tables might improve the
0070   // distributions, but also lead to more cache pollution.
0071   //
0072   // x = {3.71308, 3.44261, 3.22308, ..., 0}
0073   // f = {0.00101, 0.00266, 0.00554, ..., 1}
0074   struct Tables {
0075     double x[kMask + 2];
0076     double f[kMask + 2];
0077   };
0078   static const Tables zg_;
0079   random_internal::FastUniformBits<uint64_t> fast_u64_;
0080 };
0081
0082 }  // namespace random_internal
0083
0084 // absl::gaussian_distribution:
0085 // Generates a number conforming to a Gaussian distribution.
0086 template <typename RealType = double>
0087 class gaussian_distribution : random_internal::gaussian_distribution_base {
0088  public:
0089   using result_type = RealType;
0090
0091   class param_type {
0092    public:
0093     using distribution_type = gaussian_distribution;
0094
0095     explicit param_type(result_type mean = 0, result_type stddev = 1)
0096         : mean_(mean), stddev_(stddev) {}
0097
0098     // Returns the mean distribution parameter.  The mean specifies the location
0099     // of the peak.  The default value is 0.0.
0100     result_type mean() const { return mean_; }
0101
0102     // Returns the deviation distribution parameter.  The default value is 1.0.
0103     result_type stddev() const { return stddev_; }
0104
0105     friend bool operator==(const param_type& a, const param_type& b) {
0106       return a.mean_ == b.mean_ && a.stddev_ == b.stddev_;
0107     }
0108
0109     friend bool operator!=(const param_type& a, const param_type& b) {
0110       return !(a == b);
0111     }
0112
0113    private:
0114     result_type mean_;
0115     result_type stddev_;
0116
0117     static_assert(
0118         std::is_floating_point<RealType>::value,
0119         "Class-template absl::gaussian_distribution<> must be parameterized "
0120         "using a floating-point type.");
0121   };
0122
0123   gaussian_distribution() : gaussian_distribution(0) {}
0124
0125   explicit gaussian_distribution(result_type mean, result_type stddev = 1)
0126       : param_(mean, stddev) {}
0127
0128   explicit gaussian_distribution(const param_type& p) : param_(p) {}
0129
0130   void reset() {}
0131
0132   // Generating functions
0133   template <typename URBG>
0134   result_type operator()(URBG& g) {  // NOLINT(runtime/references)
0135     return (*this)(g, param_);
0136   }
0137
0138   template <typename URBG>
0139   result_type operator()(URBG& g,  // NOLINT(runtime/references)
0140                          const param_type& p);
0141
0142   param_type param() const { return param_; }
0143   void param(const param_type& p) { param_ = p; }
0144
0145   result_type(min)() const {
0146     return -std::numeric_limits<result_type>::infinity();
0147   }
0148   result_type(max)() const {
0149     return std::numeric_limits<result_type>::infinity();
0150   }
0151
0152   result_type mean() const { return param_.mean(); }
0153   result_type stddev() const { return param_.stddev(); }
0154
0155   friend bool operator==(const gaussian_distribution& a,
0156                          const gaussian_distribution& b) {
0157     return a.param_ == b.param_;
0158   }
0159   friend bool operator!=(const gaussian_distribution& a,
0160                          const gaussian_distribution& b) {
0161     return a.param_ != b.param_;
0162   }
0163
0164  private:
0165   param_type param_;
0166 };
0167
0168 // --------------------------------------------------------------------------
0169 // Implementation details only below
0170 // --------------------------------------------------------------------------
0171
0172 template <typename RealType>
0173 template <typename URBG>
0174 typename gaussian_distribution<RealType>::result_type
0175 gaussian_distribution<RealType>::operator()(
0176     URBG& g,  // NOLINT(runtime/references)
0177     const param_type& p) {
0178   return p.mean() + p.stddev() * static_cast<result_type>(zignor(g));
0179 }
0180
0181 template <typename CharT, typename Traits, typename RealType>
0182 std::basic_ostream<CharT, Traits>& operator<<(
0183     std::basic_ostream<CharT, Traits>& os,  // NOLINT(runtime/references)
0184     const gaussian_distribution<RealType>& x) {
0185   auto saver = random_internal::make_ostream_state_saver(os);
0186   os.precision(random_internal::stream_precision_helper<RealType>::kPrecision);
0187   os << x.mean() << os.fill() << x.stddev();
0188   return os;
0189 }
0190
0191 template <typename CharT, typename Traits, typename RealType>
0192 std::basic_istream<CharT, Traits>& operator>>(
0193     std::basic_istream<CharT, Traits>& is,  // NOLINT(runtime/references)
0194     gaussian_distribution<RealType>& x) {   // NOLINT(runtime/references)
0195   using result_type = typename gaussian_distribution<RealType>::result_type;
0196   using param_type = typename gaussian_distribution<RealType>::param_type;
0197
0198   auto saver = random_internal::make_istream_state_saver(is);
0199   auto mean = random_internal::read_floating_point<result_type>(is);
0200   if (is.fail()) return is;
0201   auto stddev = random_internal::read_floating_point<result_type>(is);
0202   if (!is.fail()) {
0203     x.param(param_type(mean, stddev));
0204   }
0205   return is;
0206 }
0207
0208 namespace random_internal {
0209
0210 template <typename URBG>
0211 inline double gaussian_distribution_base::zignor_fallback(URBG& g, bool neg) {
0212   using random_internal::GeneratePositiveTag;
0213   using random_internal::GenerateRealFromBits;
0214
0215   // This fallback path happens approximately 0.05% of the time.
0216   double x, y;
0217   do {
0218     // kRInv = 1/r, U(0, 1)
0219     x = kRInv *
0220         std::log(GenerateRealFromBits<double, GeneratePositiveTag, false>(
0221             fast_u64_(g)));
0222     y = -std::log(
0223         GenerateRealFromBits<double, GeneratePositiveTag, false>(fast_u64_(g)));
0224   } while ((y + y) < (x * x));
0225   return neg ? (x - kR) : (kR - x);
0226 }
0227
0228 template <typename URBG>
0229 inline double gaussian_distribution_base::zignor(
0230     URBG& g) {  // NOLINT(runtime/references)
0231   using random_internal::GeneratePositiveTag;
0232   using random_internal::GenerateRealFromBits;
0233   using random_internal::GenerateSignedTag;
0234
0235   while (true) {
0236     // We use a single uint64_t to generate both a double and a strip.
0237     // These bits are unused when the generated double is > 1/2^5.
0238     // This may introduce some bias from the duplicated low bits of small
0239     // values (those smaller than 1/2^5, which all end up on the left tail).
0240     uint64_t bits = fast_u64_(g);
0241     int i = static_cast<int>(bits & kMask);  // pick a random strip
0242     double j = GenerateRealFromBits<double, GenerateSignedTag, false>(
0243         bits);  // U(-1, 1)
0244     const double x = j * zg_.x[i];
0245
0246     // Retangular box. Handles >97% of all cases.
0247     // For any given box, this handles between 75% and 99% of values.
0248     // Equivalent to U(01) < (x[i+1] / x[i]), and when i == 0, ~93.5%
0249     if (std::abs(x) < zg_.x[i + 1]) {
0250       return x;
0251     }
0252
0253     // i == 0: Base box. Sample using a ratio of uniforms.
0254     if (i == 0) {
0255       // This path happens about 0.05% of the time.
0256       return zignor_fallback(g, j < 0);
0257     }
0258
0259     // i > 0: Wedge samples using precomputed values.
0260     double v = GenerateRealFromBits<double, GeneratePositiveTag, false>(
0261         fast_u64_(g));  // U(0, 1)
0262     if ((zg_.f[i + 1] + v * (zg_.f[i] - zg_.f[i + 1])) <
0263         std::exp(-0.5 * x * x)) {
0264       return x;
0265     }
0266
0267     // The wedge was missed; reject the value and try again.
0268   }
0269 }
0270
0271 }  // namespace random_internal
0272 ABSL_NAMESPACE_END
0273 }  // namespace absl
0274
0275 #endif  // ABSL_RANDOM_GAUSSIAN_DISTRIBUTION_H_