libstdc++: Optimise GCD algorithms

The current std::gcd and std::chrono::duration::_S_gcd algorithms are
both recursive. This is potentially expensive to evaluate in constant
expressions, because each level of recursion makes a new copy of the
function to evaluate. The maximum number of steps is bounded
(proportional to the number of decimal digits in the smaller value) and
so unlikely to exceed the limit for constexpr nesting, but the memory
usage is still suboptimal. By using an iterative algorithm we avoid
that compile-time cost. Because looping in constexpr functions is not
allowed until C++14, we need to keep the recursive implementation in
duration::_S_gcd for C++11 mode.

For std::gcd we can also optimise runtime performance by using the
binary GCD algorithm.

libstdc++-v3/ChangeLog:

	* include/std/chrono (duration::_S_gcd): Use iterative algorithm
	for C++14 and later.
	* include/std/numeric (__detail::__gcd): Replace recursive
	Euclidean algorithm with iterative version of binary GCD algorithm.
	* testsuite/26_numerics/gcd/1.cc: Test additional inputs.
	* testsuite/26_numerics/gcd/gcd_neg.cc: Adjust dg-error lines.
	* testsuite/26_numerics/lcm/lcm_neg.cc: Likewise.
	* testsuite/experimental/numeric/gcd.cc: Test additional inputs.
	* testsuite/26_numerics/gcd/2.cc: New test.

1 files changed, 30 insertions, 4 deletions

diff --git a/libstdc++-v3/include/std/numeric b/libstdc++-v3/include/std/numeric
index bd70a52..2de6aaf 100644
--- a/libstdc++-v3/include/std/numeric
+++ b/libstdc++-v3/include/std/numeric
@@ -76,6 +76,7 @@
 
 #if __cplusplus >= 201402L
 #include <type_traits>
+#include <bit>
 
 namespace std _GLIBCXX_VISIBILITY(default)
 {
@@ -97,15 +98,40 @@ namespace __detail
 
   template<typename _Up> void __absu(bool) = delete;
 
-  // GCD implementation
+  // GCD implementation, using Stein's algorithm
   template<typename _Tp>
     constexpr _Tp
     __gcd(_Tp __m, _Tp __n)
     {
       static_assert(is_unsigned<_Tp>::value, "type must be unsigned");
-      return __m == 0 ? __n
-	: __n == 0 ? __m
-	: __detail::__gcd(__n, _Tp(__m % __n));
+
+      if (__m == 0)
+	return __n;
+      if (__n == 0)
+	return __m;
+
+      const int __i = std::__countr_zero(__m);
+      __m >>= __i;
+      const int __j = std::__countr_zero(__n);
+      __n >>= __j;
+      const int __k = __i < __j ? __i : __j; // min(i, j)
+
+      while (true)
+	{
+	  if (__m > __n)
+	    {
+	      _Tp __tmp = __m;
+	      __m = __n;
+	      __n = __tmp;
+	    }
+
+	  __n -= __m;
+
+	  if (__n == 0)
+	    return __m << __k;
+
+	  __n >>= std::__countr_zero(__n);
+	}
     }
 
   // LCM implementation