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.

From-SVN: r272616

5 files changed, 52 insertions, 17 deletions

diff --git a/libstdc++-v3/ChangeLog b/libstdc++-v3/ChangeLog
index 2acb150..35c79f9 100644
--- a/libstdc++-v3/ChangeLog
+++ b/libstdc++-v3/ChangeLog
@@ -1,5 +1,12 @@
 2019-06-24  Jonathan Wakely  <jwakely@redhat.com>
 
+	* 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.
+
 	* testsuite/18_support/headers/cfloat/values_c++17.cc: New test.
 
 2019-06-20  Jonathan Wakely  <jwakely@redhat.com>
diff --git a/libstdc++-v3/include/std/numeric b/libstdc++-v3/include/std/numeric
index fc2242f..af68446 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 87a7498..05e2143 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 4db01ae..3a0f5bb 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 9c6e411..32a966e 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();