diff --git a/include/cmath b/include/cmath index f5f62adcf..56d17bfdd 100644 --- a/include/cmath +++ b/include/cmath @@ -304,11 +304,15 @@ long double truncl(long double x); #include <__config> #include #include +#include #if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER) #pragma GCC system_header #endif +_LIBCPP_PUSH_MACROS +#include <__undef_macros> + _LIBCPP_BEGIN_NAMESPACE_STD using ::signbit; @@ -607,6 +611,38 @@ __libcpp_isfinite_or_builtin(_A1 __lcpp_x) _NOEXCEPT return isfinite(__lcpp_x); } +template ::digits > numeric_limits<_IntT>::digits), + int _Bits = (numeric_limits<_IntT>::digits - numeric_limits<_FloatT>::digits)> +_LIBCPP_INLINE_VISIBILITY +_LIBCPP_CONSTEXPR _IntT __max_representable_int_for_float() _NOEXCEPT { + static_assert(is_floating_point<_FloatT>::value, "must be a floating point type"); + static_assert(is_integral<_IntT>::value, "must be an integral type"); + static_assert(numeric_limits<_FloatT>::radix == 2, "FloatT has incorrect radix"); + static_assert(is_same<_FloatT, float>::value || is_same<_FloatT, double>::value + || is_same<_FloatT,long double>::value, "unsupported floating point type"); + return _FloatBigger ? numeric_limits<_IntT>::max() : (numeric_limits<_IntT>::max() >> _Bits << _Bits); +} + +// Convert a floating point number to the specified integral type after +// clamping to the integral types representable range. +// +// The behavior is undefined if `__r` is NaN. +template +_LIBCPP_INLINE_VISIBILITY +_IntT __clamp_to_integral(_RealT __r) _NOEXCEPT { + using _Lim = std::numeric_limits<_IntT>; + const _IntT _MaxVal = std::__max_representable_int_for_float<_IntT, _RealT>(); + if (__r >= ::nextafter(static_cast<_RealT>(_MaxVal), INFINITY)) { + return _Lim::max(); + } else if (__r <= _Lim::lowest()) { + return _Lim::min(); + } + return static_cast<_IntT>(__r); +} + _LIBCPP_END_NAMESPACE_STD +_LIBCPP_POP_MACROS + #endif // _LIBCPP_CMATH diff --git a/include/random b/include/random index 724bd0fc2..9b29e14a6 100644 --- a/include/random +++ b/include/random @@ -4593,7 +4593,10 @@ public: template poisson_distribution<_IntType>::param_type::param_type(double __mean) - : __mean_(__mean) + // According to the standard `inf` is a valid input, but it causes the + // distribution to hang, so we replace it with the maximum representable + // mean. + : __mean_(isinf(__mean) ? numeric_limits::max() : __mean) { if (__mean_ < 10) { @@ -4611,7 +4614,7 @@ poisson_distribution<_IntType>::param_type::param_type(double __mean) { __s_ = _VSTD::sqrt(__mean_); __d_ = 6 * __mean_ * __mean_; - __l_ = static_cast(__mean_ - 1.1484); + __l_ = std::trunc(__mean_ - 1.1484); __omega_ = .3989423 / __s_; double __b1_ = .4166667E-1 / __mean_; double __b2_ = .3 * __b1_ * __b1_; @@ -4628,12 +4631,12 @@ template _IntType poisson_distribution<_IntType>::operator()(_URNG& __urng, const param_type& __pr) { - result_type __x; + double __tx; uniform_real_distribution __urd; if (__pr.__mean_ < 10) { - __x = 0; - for (double __p = __urd(__urng); __p > __pr.__l_; ++__x) + __tx = 0; + for (double __p = __urd(__urng); __p > __pr.__l_; ++__tx) __p *= __urd(__urng); } else @@ -4643,19 +4646,19 @@ poisson_distribution<_IntType>::operator()(_URNG& __urng, const param_type& __pr double __u; if (__g > 0) { - __x = static_cast(__g); - if (__x >= __pr.__l_) - return __x; - __difmuk = __pr.__mean_ - __x; + __tx = std::trunc(__g); + if (__tx >= __pr.__l_) + return std::__clamp_to_integral(__tx); + __difmuk = __pr.__mean_ - __tx; __u = __urd(__urng); if (__pr.__d_ * __u >= __difmuk * __difmuk * __difmuk) - return __x; + return std::__clamp_to_integral(__tx); } exponential_distribution __edist; for (bool __using_exp_dist = false; true; __using_exp_dist = true) { double __e; - if (__using_exp_dist || __g < 0) + if (__using_exp_dist || __g <= 0) { double __t; do @@ -4665,31 +4668,31 @@ poisson_distribution<_IntType>::operator()(_URNG& __urng, const param_type& __pr __u += __u - 1; __t = 1.8 + (__u < 0 ? -__e : __e); } while (__t <= -.6744); - __x = __pr.__mean_ + __pr.__s_ * __t; - __difmuk = __pr.__mean_ - __x; + __tx = std::trunc(__pr.__mean_ + __pr.__s_ * __t); + __difmuk = __pr.__mean_ - __tx; __using_exp_dist = true; } double __px; double __py; - if (__x < 10) + if (__tx < 10 && __tx >= 0) { const double __fac[] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880}; __px = -__pr.__mean_; - __py = _VSTD::pow(__pr.__mean_, (double)__x) / __fac[__x]; + __py = _VSTD::pow(__pr.__mean_, (double)__tx) / __fac[static_cast(__tx)]; } else { - double __del = .8333333E-1 / __x; + double __del = .8333333E-1 / __tx; __del -= 4.8 * __del * __del * __del; - double __v = __difmuk / __x; + double __v = __difmuk / __tx; if (_VSTD::abs(__v) > 0.25) - __px = __x * _VSTD::log(1 + __v) - __difmuk - __del; + __px = __tx * _VSTD::log(1 + __v) - __difmuk - __del; else - __px = __x * __v * __v * (((((((.1250060 * __v + -.1384794) * + __px = __tx * __v * __v * (((((((.1250060 * __v + -.1384794) * __v + .1421878) * __v + -.1661269) * __v + .2000118) * __v + -.2500068) * __v + .3333333) * __v + -.5) - __del; - __py = .3989423 / _VSTD::sqrt(__x); + __py = .3989423 / _VSTD::sqrt(__tx); } double __r = (0.5 - __difmuk) / __pr.__s_; double __r2 = __r * __r; @@ -4709,7 +4712,7 @@ poisson_distribution<_IntType>::operator()(_URNG& __urng, const param_type& __pr } } } - return __x; + return std::__clamp_to_integral(__tx); } template diff --git a/test/libcxx/numerics/c.math/undef_min_max.pass.cpp b/test/libcxx/numerics/c.math/undef_min_max.pass.cpp new file mode 100644 index 000000000..ba1696eb2 --- /dev/null +++ b/test/libcxx/numerics/c.math/undef_min_max.pass.cpp @@ -0,0 +1,19 @@ +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#if defined(__GNUC__) || defined(__clang__) +#pragma GCC diagnostic ignored "-W#warnings" +#endif + +#define min THIS IS A NASTY MACRO! +#define max THIS IS A NASTY MACRO! + +#include + +#include "test_macros.h" + +int main(int, char**) { return 0; } diff --git a/test/libcxx/numerics/clamp_to_integral.pass.cpp b/test/libcxx/numerics/clamp_to_integral.pass.cpp new file mode 100644 index 000000000..cb3336f52 --- /dev/null +++ b/test/libcxx/numerics/clamp_to_integral.pass.cpp @@ -0,0 +1,90 @@ +//===----------------------------------------------------------------------===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +// __clamp_to_integral(RealT) + +// Test the conversion function that truncates floating point types to the +// closest representable value for the specified integer type, or +// numeric_limits::max()/min() if the value isn't representable. + +#include +#include +#include + +template +void test() { + typedef std::numeric_limits Lim; + const bool MaxIsRepresentable = sizeof(IntT) < 8; + const bool IsSigned = std::is_signed::value; + struct TestCase { + double Input; + IntT Expect; + bool IsRepresentable; + } TestCases[] = { + {0, 0, true}, + {1, 1, true}, + {IsSigned ? static_cast(-1) : 0, + IsSigned ? static_cast(-1) : 0, true}, + {Lim::lowest(), Lim::lowest(), true}, + {static_cast(Lim::max()), Lim::max(), MaxIsRepresentable}, + {static_cast(Lim::max()) + 1, Lim::max(), false}, + {static_cast(Lim::max()) + 1024, Lim::max(), false}, + {nextafter(static_cast(Lim::max()), INFINITY), Lim::max(), false}, + }; + for (TestCase TC : TestCases) { + auto res = std::__clamp_to_integral(TC.Input); + assert(res == TC.Expect); + if (TC.IsRepresentable) { + auto other = static_cast(std::trunc(TC.Input)); + assert(res == other); + } else + assert(res == Lim::min() || res == Lim::max()); + } +} + +template +void test_float() { + typedef std::numeric_limits Lim; + const bool MaxIsRepresentable = sizeof(IntT) < 4; + ((void)MaxIsRepresentable); + const bool IsSigned = std::is_signed::value; + struct TestCase { + float Input; + IntT Expect; + bool IsRepresentable; + } TestCases[] = { + {0, 0, true}, + {1, 1, true}, + {IsSigned ? static_cast(-1) : 0, + IsSigned ? static_cast(-1) : 0, true}, + {Lim::lowest(), Lim::lowest(), true}, + {static_cast(Lim::max()), Lim::max(), MaxIsRepresentable }, + {nextafter(static_cast(Lim::max()), INFINITY), Lim::max(), false}, + }; + for (TestCase TC : TestCases) { + auto res = std::__clamp_to_integral(TC.Input); + assert(res == TC.Expect); + if (TC.IsRepresentable) { + auto other = static_cast(std::trunc(TC.Input)); + assert(res == other); + } else + assert(res == Lim::min() || res == Lim::max()); + } +} + +int main() { + test(); + test(); + test(); + test(); + test(); + test(); + test_float(); + test_float(); + test_float(); +} diff --git a/test/std/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.geo/eval.pass.cpp b/test/std/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.geo/eval.pass.cpp index 6e6072a4f..3c726a1db 100644 --- a/test/std/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.geo/eval.pass.cpp +++ b/test/std/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.geo/eval.pass.cpp @@ -29,6 +29,20 @@ sqr(T x) return x * x; } +struct Eng : std::mt19937 { + using Base = std::mt19937; + using Base::Base; +}; + +void test_small_inputs() { + Eng engine; + std::geometric_distribution distribution(5.45361e-311); + for (auto i=0; i < 1000; ++i) { + volatile auto res = distribution(engine); + ((void)res); + } +} + void test1() { @@ -295,4 +309,5 @@ int main() test4(); test5(); test6(); + test_small_inputs(); } diff --git a/test/std/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.poisson/eval.pass.cpp b/test/std/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.poisson/eval.pass.cpp index 12fcfa354..e7127d773 100644 --- a/test/std/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.poisson/eval.pass.cpp +++ b/test/std/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.poisson/eval.pass.cpp @@ -29,6 +29,68 @@ sqr(T x) return x * x; } +void test_bad_ranges() { + // Test cases where the mean is around the largest representable integer for + // `result_type`. These cases don't generate valid poisson distributions, but + // at least they don't blow up. + std::mt19937 eng; + + { + std::poisson_distribution distribution(32710.9); + for (int i=0; i < 1000; ++i) { + volatile std::int16_t res = distribution(eng); + ((void)res); + } + } + { + std::poisson_distribution distribution(std::numeric_limits::max()); + for (int i=0; i < 1000; ++i) { + volatile std::int16_t res = distribution(eng); + ((void)res); + } + } + { + std::poisson_distribution distribution( + static_cast(std::numeric_limits::max()) + 10); + for (int i=0; i < 1000; ++i) { + volatile std::int16_t res = distribution(eng); + ((void)res); + } + } + { + std::poisson_distribution distribution( + static_cast(std::numeric_limits::max()) * 2); + for (int i=0; i < 1000; ++i) { + volatile std::int16_t res = distribution(eng); + ((void)res); + } + } + { + // We convert `INF` to `DBL_MAX` otherwise the distribution will hang. + std::poisson_distribution distribution(std::numeric_limits::infinity()); + for (int i=0; i < 1000; ++i) { + volatile std::int16_t res = distribution(eng); + ((void)res); + } + } + { + std::poisson_distribution distribution(0); + for (int i=0; i < 1000; ++i) { + volatile std::int16_t res = distribution(eng); + ((void)res); + } + } + { + // We convert `INF` to `DBL_MAX` otherwise the distribution will hang. + std::poisson_distribution distribution(-100); + for (int i=0; i < 1000; ++i) { + volatile std::int16_t res = distribution(eng); + ((void)res); + } + } +} + + int main() { { @@ -148,4 +210,6 @@ int main() assert(std::abs((skew - x_skew) / x_skew) < 0.01); assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01); } + + test_bad_ranges(); }