From 8e5d0aa8b87cebe5066578a76ec444e2554fa41f Mon Sep 17 00:00:00 2001 From: Igor Zhukov Date: Mon, 20 Jun 2022 07:46:08 +0700 Subject: [PATCH] Recalculate lerp if we got infinity. Eliminates some overflows. (#1918) Co-authored-by: Adam Bucior <35536269+AdamBucior@users.noreply.github.com> Co-authored-by: Alexander Bolz Co-authored-by: Curtis J Bezault Co-authored-by: Matt Stephanson Co-authored-by: Michael Schellenberger Costa Co-authored-by: statementreply Co-authored-by: Stephan T. Lavavej --- stl/inc/cmath | 27 ++++++- .../std/tests/P0811R3_midpoint_lerp/test.cpp | 81 +++++++++++++++++++ 2 files changed, 107 insertions(+), 1 deletion(-) diff --git a/stl/inc/cmath b/stl/inc/cmath index c1158d26d..21d8fe1d9 100644 --- a/stl/inc/cmath +++ b/stl/inc/cmath @@ -1335,6 +1335,31 @@ _NODISCARD auto hypot(const _Ty1 _Dx, const _Ty2 _Dy, const _Ty3 _Dz) { } #if _HAS_CXX20 +template +_NODISCARD constexpr _Ty _Linear_for_lerp(const _Ty _ArgA, const _Ty _ArgB, const _Ty _ArgT) { + if (_STD is_constant_evaluated()) { + auto _Smaller = _ArgT; + auto _Larger = _ArgB - _ArgA; + auto _Abs_smaller = _Float_abs(_Smaller); + auto _Abs_larger = _Float_abs(_Larger); + if (_Abs_larger < _Abs_smaller) { + _STD swap(_Smaller, _Larger); + _STD swap(_Abs_smaller, _Abs_larger); + } + + if (_Abs_smaller > 1) { + // _Larger is too large to be subnormal, so scaling by 0.5 is exact, and the product _Smaller * _Larger is + // large enough that if _ArgA is subnormal, it will be too small to contribute anyway and this way can + // sometimes avoid overflow problems. + return 2 * (_Ty{0.5} * _ArgA + _Smaller * (_Ty{0.5} * _Larger)); + } else { + return _ArgA + _Smaller * _Larger; + } + } + + return _STD fma(_ArgT, _ArgB - _ArgA, _ArgA); +} + template _NODISCARD constexpr _Ty _Common_lerp(const _Ty _ArgA, const _Ty _ArgB, const _Ty _ArgT) noexcept { // on a line intersecting {(0.0, _ArgA), (1.0, _ArgB)}, return the Y value for X == _ArgT @@ -1353,7 +1378,7 @@ _NODISCARD constexpr _Ty _Common_lerp(const _Ty _ArgA, const _Ty _ArgB, const _T } // exact at _ArgT == 0, monotonic except near _ArgT == 1, bounded, determinate, and consistent: - const auto _Candidate = _ArgA + _ArgT * (_ArgB - _ArgA); + const auto _Candidate = _Linear_for_lerp(_ArgA, _ArgB, _ArgT); // monotonic near _ArgT == 1: if ((_ArgT > 1) == (_ArgB > _ArgA)) { if (_ArgB > _Candidate) { diff --git a/tests/std/tests/P0811R3_midpoint_lerp/test.cpp b/tests/std/tests/P0811R3_midpoint_lerp/test.cpp index 66756569d..da12486f6 100644 --- a/tests/std/tests/P0811R3_midpoint_lerp/test.cpp +++ b/tests/std/tests/P0811R3_midpoint_lerp/test.cpp @@ -1023,6 +1023,86 @@ bool test_lerp() { return true; } +void test_gh_1917() { + // GH-1917 : lerp(1e+308, 5e+307, 4.0) spuriously overflows + using bit_type = unsigned long long; + using float_bit_type = unsigned int; + STATIC_ASSERT(bit_cast(lerp(1e+308, 5e+307, 4.0)) == bit_cast(-1e+308)); + { + ExceptGuard except; + + assert(bit_cast(lerp(1e+308, 5e+307, 4.0)) == bit_cast(-1e+308)); + assert(check_feexcept(0)); + } + STATIC_ASSERT(bit_cast(lerp(2e+38f, 1e+38f, 4.0f)) == bit_cast(-2e+38f)); + { + ExceptGuard except; + + assert(bit_cast(lerp(2e+38f, 1e+38f, 4.0f)) == bit_cast(-2e+38f)); + assert(check_feexcept(0)); + } +#ifdef _M_FP_STRICT + { + ExceptGuard except; + RoundGuard round{FE_UPWARD}; + + assert(bit_cast(lerp(1e+308, 5e+307, 4.0)) == bit_cast(-1e+308)); + assert(check_feexcept(0)); + } + { + ExceptGuard except; + RoundGuard round{FE_UPWARD}; + + assert(bit_cast(lerp(2e+38f, 1e+38f, 4.0f)) == bit_cast(-2e+38f)); + assert(check_feexcept(0)); + } + { + ExceptGuard except; + RoundGuard round{FE_DOWNWARD}; + + assert(bit_cast(lerp(1e+308, 5e+307, 4.0)) == bit_cast(-1e+308)); + assert(check_feexcept(0)); + } + { + ExceptGuard except; + RoundGuard round{FE_DOWNWARD}; + + assert(bit_cast(lerp(2e+38f, 1e+38f, 4.0f)) == bit_cast(-2e+38f)); + assert(check_feexcept(0)); + } + { + ExceptGuard except; + RoundGuard round{FE_TOWARDZERO}; + + assert(bit_cast(lerp(1e+308, 5e+307, 4.0)) == bit_cast(-1e+308)); + assert(check_feexcept(0)); + } + { + ExceptGuard except; + RoundGuard round{FE_TOWARDZERO}; + + assert(bit_cast(lerp(2e+38f, 1e+38f, 4.0f)) == bit_cast(-2e+38f)); + assert(check_feexcept(0)); + } + { + ExceptGuard except; + const int r = feraiseexcept(FE_OVERFLOW); + + assert(r == 0); + assert(bit_cast(lerp(1e+308, 5e+307, 4.0)) == bit_cast(-1e+308)); + assert(check_feexcept(FE_OVERFLOW)); + } + { + ExceptGuard except; + const int r = feraiseexcept(FE_OVERFLOW); + + assert(r == 0); + assert(bit_cast(lerp(2e+38f, 1e+38f, 4.0f)) == bit_cast(-2e+38f)); + assert(check_feexcept(FE_OVERFLOW)); + } +#endif // _M_FP_STRICT +} + constexpr bool test_gh_2112() { // GH-2112 : std::lerp is missing Arithmetic overloads assert(lerp(0, 0, 0) == 0.0); @@ -1108,6 +1188,7 @@ int main() { test_lerp(); test_lerp(); + test_gh_1917(); test_gh_2112(); STATIC_ASSERT(test_gh_2112()); }