On 12/03/19 23:04 +0000, Jonathan Wakely wrote:
On 12/03/19 22:49 +0000, Joseph Myers wrote:
On Tue, 5 Mar 2019, Jonathan Wakely wrote:
The midpoint and lerp functions for floating point types come straight
from the P0811R3 proposal, with no attempt at optimization.
I don't know whether P0811R3 states different requirements from the public
P0811R2, but the implementation of midpoint using isnormal does *not*
satisfy "at most one inexact operation occurs" and is *not* correctly
rounded, contrary to the claims made in P0811R2.
I did wonder how the implementation in the paper was meant to meet the
stated requirements, but I didn't wonder too hard.
Consider e.g. midpoint(DBL_MIN + DBL_TRUE_MIN, DBL_MIN + DBL_TRUE_MIN).
The value DBL_MIN + DBL_TRUE_MIN is normal, but dividing it by 2 is
inexact (and so that midpoint implementation would produce DBL_MIN as
result, so failing to satisfy midpoint(x, x) == x).
Replacing isnormal(x) by something like isgreaterequal(fabs(x), MIN*2)
would avoid those inexact divisions, but there would still be spurious
overflows in non-default rounding modes for e.g. midpoint(DBL_MAX,
DBL_TRUE_MIN) in FE_UPWARD mode, so failing "No overflow occurs" if that's
meant to apply in all rounding modes.
Thanks for this review, and the useful cases to test. Ed is working on
adding some more tests, so maybe he can also look at improving the
code :-)
I've committed r272616 to make this case work. This is the proposal
author's most recent suggestion for the implementation.
Tested x86_64-linux, committed to trunk.
commit e693fb3d93cfe938d700512e8bfe70e0a5c0dd8a
Author: redi <redi@138bc75d-0d04-0410-961f-82ee72b054a4>
Date: Mon Jun 24 12:09:51 2019 +0000
Fix std::midpoint for denormal values
* include/std/numeric (midpoint(T, T)): Change implementation for
floating-point types to avoid incorrect rounding of denormals.
* testsuite/26_numerics/midpoint/floating.cc: Add check for correct
rounding with denormals.
* testsuite/26_numerics/gcd/gcd_neg.cc: Adjust dg-error line numbers.
* testsuite/26_numerics/lcm/lcm_neg.cc: Likewise.
git-svn-id: svn+ssh://gcc.gnu.org/svn/gcc/trunk@272616 138bc75d-0d04-0410-961f-82ee72b054a4
diff --git a/libstdc++-v3/include/std/numeric b/libstdc++-v3/include/std/numeric
index fc2242f3de6..af684469769 100644
--- a/libstdc++-v3/include/std/numeric
+++ b/libstdc++-v3/include/std/numeric
@@ -69,9 +69,8 @@
* @defgroup numerics Numerics
*
* Components for performing numeric operations. Includes support for
- * for complex number types, random number generation, numeric
- * (n-at-a-time) arrays, generalized numeric algorithms, and special
- * math functions.
+ * complex number types, random number generation, numeric (n-at-a-time)
+ * arrays, generalized numeric algorithms, and mathematical special functions.
*/
#if __cplusplus >= 201402L
@@ -156,11 +155,22 @@ namespace __detail
#endif // C++17
+_GLIBCXX_END_NAMESPACE_VERSION
+} // namespace std
+
+#endif // C++14
+
#if __cplusplus > 201703L
+#include <limits>
+#include <bits/std_abs.h>
+
+namespace std _GLIBCXX_VISIBILITY(default)
+{
+_GLIBCXX_BEGIN_NAMESPACE_VERSION
// midpoint
# define __cpp_lib_interpolate 201902L
-template<typename _Tp>
+ template<typename _Tp>
constexpr
enable_if_t<__and_v<is_arithmetic<_Tp>, is_same<remove_cv_t<_Tp>, _Tp>,
__not_<is_same<_Tp, bool>>>,
@@ -182,11 +192,17 @@ template<typename _Tp>
}
return __a + __k * _Tp(_Up(__M - __m) / 2);
}
- else
+ else // is_floating
{
- return __builtin_isnormal(__a) && __builtin_isnormal(__b)
- ? __a / 2 + __b / 2
- : (__a + __b) / 2;
+ constexpr _Tp __lo = numeric_limits<_Tp>::min() * 2;
+ constexpr _Tp __hi = numeric_limits<_Tp>::max() / 2;
+ if (std::abs(__a) <= __hi && std::abs(__b) <= __hi) [[likely]]
+ return (__a + __b) / 2; // always correctly rounded
+ if (std::abs(__a) < __lo) // not safe to halve __a
+ return __a + __b/2;
+ if (std::abs(__b) < __lo) // not safe to halve __b
+ return __a/2 + __b;
+ return __a/2 + __b/2; // otherwise correctly rounded
}
}
@@ -197,12 +213,10 @@ template<typename _Tp>
{
return __a + (__b - __a) / 2;
}
-#endif // C++20
-
_GLIBCXX_END_NAMESPACE_VERSION
} // namespace std
-#endif // C++14
+#endif // C++20
#if __cplusplus > 201402L
#include <bits/stl_function.h>
diff --git a/libstdc++-v3/testsuite/26_numerics/gcd/gcd_neg.cc b/libstdc++-v3/testsuite/26_numerics/gcd/gcd_neg.cc
index 87a74988fa4..05e21431313 100644
--- a/libstdc++-v3/testsuite/26_numerics/gcd/gcd_neg.cc
+++ b/libstdc++-v3/testsuite/26_numerics/gcd/gcd_neg.cc
@@ -46,9 +46,9 @@ test01()
std::gcd<const int&, const int&>(0.1, 0.1); // { dg-error "from here" }
}
+// { dg-error "integers" "" { target *-*-* } 133 }
// { dg-error "integers" "" { target *-*-* } 134 }
-// { dg-error "integers" "" { target *-*-* } 135 }
-// { dg-error "not bools" "" { target *-*-* } 136 }
-// { dg-error "not bools" "" { target *-*-* } 138 }
+// { dg-error "not bools" "" { target *-*-* } 135 }
+// { dg-error "not bools" "" { target *-*-* } 137 }
// { dg-prune-output "deleted function" }
// { dg-prune-output "invalid operands" }
diff --git a/libstdc++-v3/testsuite/26_numerics/lcm/lcm_neg.cc b/libstdc++-v3/testsuite/26_numerics/lcm/lcm_neg.cc
index 4db01ae3647..3a0f5bbe3eb 100644
--- a/libstdc++-v3/testsuite/26_numerics/lcm/lcm_neg.cc
+++ b/libstdc++-v3/testsuite/26_numerics/lcm/lcm_neg.cc
@@ -46,9 +46,9 @@ test01()
std::lcm<const int&, const int&>(0.1, 0.1); // { dg-error "from here" }
}
+// { dg-error "integers" "" { target *-*-* } 147 }
// { dg-error "integers" "" { target *-*-* } 148 }
-// { dg-error "integers" "" { target *-*-* } 149 }
-// { dg-error "not bools" "" { target *-*-* } 150 }
-// { dg-error "not bools" "" { target *-*-* } 152 }
+// { dg-error "not bools" "" { target *-*-* } 149 }
+// { dg-error "not bools" "" { target *-*-* } 151 }
// { dg-prune-output "deleted function" }
// { dg-prune-output "invalid operands" }
diff --git a/libstdc++-v3/testsuite/26_numerics/midpoint/floating.cc b/libstdc++-v3/testsuite/26_numerics/midpoint/floating.cc
index 9c6e4113b18..32a966e2c63 100644
--- a/libstdc++-v3/testsuite/26_numerics/midpoint/floating.cc
+++ b/libstdc++-v3/testsuite/26_numerics/midpoint/floating.cc
@@ -20,6 +20,7 @@
#include <numeric>
#include <limits>
+#include <cfloat>
#include <testsuite_hooks.h>
void
@@ -57,6 +58,19 @@ test03()
VERIFY( std::midpoint(9e9l, -9e9l) == 0.0l );
}
+namespace test04
+{
+ // https://gcc.gnu.org/ml/libstdc++/2019-03/msg00065.html
+ constexpr double d = DBL_MIN + DBL_TRUE_MIN;
+ static_assert( std::midpoint(d, d) == d );
+
+ constexpr float f = FLT_MIN + FLT_TRUE_MIN;
+ static_assert( std::midpoint(f, f) == f );
+
+ constexpr long double l = LDBL_MIN + LDBL_TRUE_MIN;
+ static_assert( std::midpoint(l, l) == l );
+}
+
int main()
{
test01();
