exponential_distribution.h 5.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165
  1. // Copyright 2017 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. #ifndef Y_ABSL_RANDOM_EXPONENTIAL_DISTRIBUTION_H_
  15. #define Y_ABSL_RANDOM_EXPONENTIAL_DISTRIBUTION_H_
  16. #include <cassert>
  17. #include <cmath>
  18. #include <istream>
  19. #include <limits>
  20. #include <type_traits>
  21. #include "y_absl/meta/type_traits.h"
  22. #include "y_absl/random/internal/fast_uniform_bits.h"
  23. #include "y_absl/random/internal/generate_real.h"
  24. #include "y_absl/random/internal/iostream_state_saver.h"
  25. namespace y_absl {
  26. Y_ABSL_NAMESPACE_BEGIN
  27. // y_absl::exponential_distribution:
  28. // Generates a number conforming to an exponential distribution and is
  29. // equivalent to the standard [rand.dist.pois.exp] distribution.
  30. template <typename RealType = double>
  31. class exponential_distribution {
  32. public:
  33. using result_type = RealType;
  34. class param_type {
  35. public:
  36. using distribution_type = exponential_distribution;
  37. explicit param_type(result_type lambda = 1) : lambda_(lambda) {
  38. assert(lambda > 0);
  39. neg_inv_lambda_ = -result_type(1) / lambda_;
  40. }
  41. result_type lambda() const { return lambda_; }
  42. friend bool operator==(const param_type& a, const param_type& b) {
  43. return a.lambda_ == b.lambda_;
  44. }
  45. friend bool operator!=(const param_type& a, const param_type& b) {
  46. return !(a == b);
  47. }
  48. private:
  49. friend class exponential_distribution;
  50. result_type lambda_;
  51. result_type neg_inv_lambda_;
  52. static_assert(
  53. std::is_floating_point<RealType>::value,
  54. "Class-template y_absl::exponential_distribution<> must be parameterized "
  55. "using a floating-point type.");
  56. };
  57. exponential_distribution() : exponential_distribution(1) {}
  58. explicit exponential_distribution(result_type lambda) : param_(lambda) {}
  59. explicit exponential_distribution(const param_type& p) : param_(p) {}
  60. void reset() {}
  61. // Generating functions
  62. template <typename URBG>
  63. result_type operator()(URBG& g) { // NOLINT(runtime/references)
  64. return (*this)(g, param_);
  65. }
  66. template <typename URBG>
  67. result_type operator()(URBG& g, // NOLINT(runtime/references)
  68. const param_type& p);
  69. param_type param() const { return param_; }
  70. void param(const param_type& p) { param_ = p; }
  71. result_type(min)() const { return 0; }
  72. result_type(max)() const {
  73. return std::numeric_limits<result_type>::infinity();
  74. }
  75. result_type lambda() const { return param_.lambda(); }
  76. friend bool operator==(const exponential_distribution& a,
  77. const exponential_distribution& b) {
  78. return a.param_ == b.param_;
  79. }
  80. friend bool operator!=(const exponential_distribution& a,
  81. const exponential_distribution& b) {
  82. return a.param_ != b.param_;
  83. }
  84. private:
  85. param_type param_;
  86. random_internal::FastUniformBits<uint64_t> fast_u64_;
  87. };
  88. // --------------------------------------------------------------------------
  89. // Implementation details follow
  90. // --------------------------------------------------------------------------
  91. template <typename RealType>
  92. template <typename URBG>
  93. typename exponential_distribution<RealType>::result_type
  94. exponential_distribution<RealType>::operator()(
  95. URBG& g, // NOLINT(runtime/references)
  96. const param_type& p) {
  97. using random_internal::GenerateNegativeTag;
  98. using random_internal::GenerateRealFromBits;
  99. using real_type =
  100. y_absl::conditional_t<std::is_same<RealType, float>::value, float, double>;
  101. const result_type u = GenerateRealFromBits<real_type, GenerateNegativeTag,
  102. false>(fast_u64_(g)); // U(-1, 0)
  103. // log1p(-x) is mathematically equivalent to log(1 - x) but has more
  104. // accuracy for x near zero.
  105. return p.neg_inv_lambda_ * std::log1p(u);
  106. }
  107. template <typename CharT, typename Traits, typename RealType>
  108. std::basic_ostream<CharT, Traits>& operator<<(
  109. std::basic_ostream<CharT, Traits>& os, // NOLINT(runtime/references)
  110. const exponential_distribution<RealType>& x) {
  111. auto saver = random_internal::make_ostream_state_saver(os);
  112. os.precision(random_internal::stream_precision_helper<RealType>::kPrecision);
  113. os << x.lambda();
  114. return os;
  115. }
  116. template <typename CharT, typename Traits, typename RealType>
  117. std::basic_istream<CharT, Traits>& operator>>(
  118. std::basic_istream<CharT, Traits>& is, // NOLINT(runtime/references)
  119. exponential_distribution<RealType>& x) { // NOLINT(runtime/references)
  120. using result_type = typename exponential_distribution<RealType>::result_type;
  121. using param_type = typename exponential_distribution<RealType>::param_type;
  122. result_type lambda;
  123. auto saver = random_internal::make_istream_state_saver(is);
  124. lambda = random_internal::read_floating_point<result_type>(is);
  125. if (!is.fail()) {
  126. x.param(param_type(lambda));
  127. }
  128. return is;
  129. }
  130. Y_ABSL_NAMESPACE_END
  131. } // namespace y_absl
  132. #endif // Y_ABSL_RANDOM_EXPONENTIAL_DISTRIBUTION_H_