Recalculate lerp if we got infinity. Eliminates some overflows. (microsoft#1918)

Co-authored-by: Adam Bucior <[email protected]>
Co-authored-by: Alexander Bolz <[email protected]>
Co-authored-by: Curtis J Bezault <[email protected]>
Co-authored-by: Matt Stephanson <[email protected]>
Co-authored-by: Michael Schellenberger Costa <[email protected]>
Co-authored-by: statementreply <[email protected]>
Co-authored-by: Stephan T. Lavavej <[email protected]>

Author: 8 people

stl/inc/cmath

@@ -1335,6 +1335,31 @@ _NODISCARD auto hypot(const _Ty1 _Dx, const _Ty2 _Dy, const _Ty3 _Dz) {
 }
 
 #if _HAS_CXX20
+template <class _Ty>
+_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 <class _Ty>
 _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) {

tests/std/tests/P0811R3_midpoint_lerp/test.cpp

@@ -1023,6 +1023,86 @@ bool test_lerp() {
     return true;
 }
 
+void test_gh_1917() {
+    // GH-1917 <cmath>: 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<bit_type>(lerp(1e+308, 5e+307, 4.0)) == bit_cast<bit_type>(-1e+308));
+    {
+        ExceptGuard except;
+
+        assert(bit_cast<bit_type>(lerp(1e+308, 5e+307, 4.0)) == bit_cast<bit_type>(-1e+308));
+        assert(check_feexcept(0));
+    }
+    STATIC_ASSERT(bit_cast<float_bit_type>(lerp(2e+38f, 1e+38f, 4.0f)) == bit_cast<float_bit_type>(-2e+38f));
+    {
+        ExceptGuard except;
+
+        assert(bit_cast<float_bit_type>(lerp(2e+38f, 1e+38f, 4.0f)) == bit_cast<float_bit_type>(-2e+38f));
+        assert(check_feexcept(0));
+    }
+#ifdef _M_FP_STRICT
+    {
+        ExceptGuard except;
+        RoundGuard round{FE_UPWARD};
+
+        assert(bit_cast<bit_type>(lerp(1e+308, 5e+307, 4.0)) == bit_cast<bit_type>(-1e+308));
+        assert(check_feexcept(0));
+    }
+    {
+        ExceptGuard except;
+        RoundGuard round{FE_UPWARD};
+
+        assert(bit_cast<float_bit_type>(lerp(2e+38f, 1e+38f, 4.0f)) == bit_cast<float_bit_type>(-2e+38f));
+        assert(check_feexcept(0));
+    }
+    {
+        ExceptGuard except;
+        RoundGuard round{FE_DOWNWARD};
+
+        assert(bit_cast<bit_type>(lerp(1e+308, 5e+307, 4.0)) == bit_cast<bit_type>(-1e+308));
+        assert(check_feexcept(0));
+    }
+    {
+        ExceptGuard except;
+        RoundGuard round{FE_DOWNWARD};
+
+        assert(bit_cast<float_bit_type>(lerp(2e+38f, 1e+38f, 4.0f)) == bit_cast<float_bit_type>(-2e+38f));
+        assert(check_feexcept(0));
+    }
+    {
+        ExceptGuard except;
+        RoundGuard round{FE_TOWARDZERO};
+
+        assert(bit_cast<bit_type>(lerp(1e+308, 5e+307, 4.0)) == bit_cast<bit_type>(-1e+308));
+        assert(check_feexcept(0));
+    }
+    {
+        ExceptGuard except;
+        RoundGuard round{FE_TOWARDZERO};
+
+        assert(bit_cast<float_bit_type>(lerp(2e+38f, 1e+38f, 4.0f)) == bit_cast<float_bit_type>(-2e+38f));
+        assert(check_feexcept(0));
+    }
+    {
+        ExceptGuard except;
+        const int r = feraiseexcept(FE_OVERFLOW);
+
+        assert(r == 0);
+        assert(bit_cast<bit_type>(lerp(1e+308, 5e+307, 4.0)) == bit_cast<bit_type>(-1e+308));
+        assert(check_feexcept(FE_OVERFLOW));
+    }
+    {
+        ExceptGuard except;
+        const int r = feraiseexcept(FE_OVERFLOW);
+
+        assert(r == 0);
+        assert(bit_cast<float_bit_type>(lerp(2e+38f, 1e+38f, 4.0f)) == bit_cast<float_bit_type>(-2e+38f));
+        assert(check_feexcept(FE_OVERFLOW));
+    }
+#endif // _M_FP_STRICT
+}
+
 constexpr bool test_gh_2112() {
     // GH-2112 <cmath>: std::lerp is missing Arithmetic overloads
     assert(lerp(0, 0, 0) == 0.0);
@@ -1108,6 +1188,7 @@ int main() {
     test_lerp<double>();
     test_lerp<long double>();
 
+    test_gh_1917();
     test_gh_2112();
     STATIC_ASSERT(test_gh_2112());
 }