[libc++] Implement std::gcd using the binary version

The binary version is four times faster than current implementation
in my setup, and generally considered a better implementation.

Code inspired by https://en.algorithmica.org/hpc/algorithms/gcd/
which itself is inspired by https://lemire.me/blog/2013/12/26/fastest-way-to-compute-the-greatest-common-divisor/

Fix #77648

Author: serge-sans-paille

libcxx/include/__bit/countr.h

@@ -12,6 +12,7 @@
 #include <__bit/rotate.h>
 #include <__concepts/arithmetic.h>
 #include <__config>
+#include <__type_traits/is_unsigned_integer.h>
 #include <limits>
 
 #if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER)
@@ -35,10 +36,11 @@ _LIBCPP_NODISCARD inline _LIBCPP_HIDE_FROM_ABI _LIBCPP_CONSTEXPR int __libcpp_ct
   return __builtin_ctzll(__x);
 }
 
-#if _LIBCPP_STD_VER >= 20
+#if _LIBCPP_STD_VER >= 17
 
-template <__libcpp_unsigned_integer _Tp>
-_LIBCPP_NODISCARD_EXT _LIBCPP_HIDE_FROM_ABI constexpr int countr_zero(_Tp __t) noexcept {
+template <typename _Tp>
+_LIBCPP_NODISCARD_EXT _LIBCPP_HIDE_FROM_ABI _LIBCPP_CONSTEXPR int __countr_zero(_Tp __t) _NOEXCEPT {
+  static_assert(__libcpp_is_unsigned_integer<_Tp>::value, "__countr_zero requires an unsigned integer type");
   if (__t == 0)
     return numeric_limits<_Tp>::digits;
 
@@ -59,6 +61,15 @@ _LIBCPP_NODISCARD_EXT _LIBCPP_HIDE_FROM_ABI constexpr int countr_zero(_Tp __t) n
   }
 }
 
+#endif
+
+#if _LIBCPP_STD_VER >= 20
+
+template <__libcpp_unsigned_integer _Tp>
+_LIBCPP_NODISCARD_EXT _LIBCPP_HIDE_FROM_ABI constexpr int countr_zero(_Tp __t) noexcept {
+  return std::__countr_zero(__t);
+}
+
 template <__libcpp_unsigned_integer _Tp>
 _LIBCPP_NODISCARD_EXT _LIBCPP_HIDE_FROM_ABI constexpr int countr_one(_Tp __t) noexcept {
   return __t != numeric_limits<_Tp>::max() ? std::countr_zero(static_cast<_Tp>(~__t)) : numeric_limits<_Tp>::digits;

libcxx/include/__numeric/gcd_lcm.h

@@ -10,7 +10,9 @@
 #ifndef _LIBCPP___NUMERIC_GCD_LCM_H
 #define _LIBCPP___NUMERIC_GCD_LCM_H
 
+#include <__algorithm/min.h>
 #include <__assert>
+#include <__bit/countr.h>
 #include <__config>
 #include <__type_traits/common_type.h>
 #include <__type_traits/is_integral.h>
@@ -50,9 +52,25 @@ struct __ct_abs<_Result, _Source, false> {
 };
 
 template <class _Tp>
-_LIBCPP_CONSTEXPR _LIBCPP_HIDDEN _Tp __gcd(_Tp __m, _Tp __n) {
+_LIBCPP_CONSTEXPR _LIBCPP_HIDDEN _Tp __gcd(_Tp __a, _Tp __b) {
   static_assert((!is_signed<_Tp>::value), "");
-  return __n == 0 ? __m : std::__gcd<_Tp>(__n, __m % __n);
+  if (__a == 0)
+    return __b;
+  if (__b == 0)
+    return __a;
+
+  int __az    = std::__countr_zero(__a);
+  int __bz    = std::__countr_zero(__b);
+  int __shift = std::min(__az, __bz);
+  __b >>= __bz;
+  while (__a != 0) {
+    __a >>= __az;
+    _Tp __absdiff = __a > __b ? __a - __b : __b - __a;
+    __b           = std::min(__a, __b);
+    __a           = __absdiff;
+    __az          = std::__countr_zero(__absdiff);
+  }
+  return __b << __shift;
 }
 
 template <class _Tp, class _Up>

libcxx/test/libcxx/transitive_includes/cxx03.csv

@@ -549,6 +549,7 @@ numeric cstddef
 numeric cstdint
 numeric execution
 numeric functional
+numeric initializer_list
 numeric iterator
 numeric limits
 numeric new

libcxx/test/libcxx/transitive_includes/cxx11.csv

@@ -554,6 +554,7 @@ numeric cstddef
 numeric cstdint
 numeric execution
 numeric functional
+numeric initializer_list
 numeric iterator
 numeric limits
 numeric new

libcxx/test/libcxx/transitive_includes/cxx14.csv

@@ -556,6 +556,7 @@ numeric cstddef
 numeric cstdint
 numeric execution
 numeric functional
+numeric initializer_list
 numeric iterator
 numeric limits
 numeric new

libcxx/test/libcxx/transitive_includes/cxx17.csv

@@ -556,6 +556,7 @@ numeric cstddef
 numeric cstdint
 numeric execution
 numeric functional
+numeric initializer_list
 numeric iterator
 numeric limits
 numeric new

libcxx/test/libcxx/transitive_includes/cxx20.csv

@@ -561,6 +561,7 @@ numeric cstddef
 numeric cstdint
 numeric execution
 numeric functional
+numeric initializer_list
 numeric iterator
 numeric limits
 numeric new

libcxx/test/libcxx/transitive_includes/cxx26.csv

@@ -173,6 +173,29 @@ experimental/simd limits
 experimental/type_traits initializer_list
 experimental/type_traits type_traits
 experimental/utility utility
+experimental/vector experimental/memory_resource
+experimental/vector vector
+ext/hash_map algorithm
+ext/hash_map cmath
+ext/hash_map cstddef
+ext/hash_map cstdint
+ext/hash_map cstring
+ext/hash_map functional
+ext/hash_map initializer_list
+ext/hash_map limits
+ext/hash_map new
+ext/hash_map stdexcept
+ext/hash_map string
+ext/hash_set algorithm
+ext/hash_set cmath
+ext/hash_set cstddef
+ext/hash_set cstdint
+ext/hash_set cstring
+ext/hash_set functional
+ext/hash_set initializer_list
+ext/hash_set limits
+ext/hash_set new
+ext/hash_set string
 filesystem compare
 filesystem cstddef
 filesystem cstdint

libcxx/test/std/numerics/numeric.ops/numeric.ops.gcd/gcd.pass.cpp

@@ -17,6 +17,7 @@
 #include <cassert>
 #include <climits>
 #include <cstdint>
+#include <random>
 #include <type_traits>
 
 #include "test_macros.h"
@@ -48,6 +49,68 @@ constexpr bool test0(int in1, int in2, int out)
     return true;
 }
 
+template <typename T>
+T basic_gcd_(T m, T n) {
+  return n == 0 ? m : basic_gcd_<T>(n, m % n);
+}
+
+template <typename T>
+T basic_gcd(T m, T n) {
+  using Tp = std::make_unsigned_t<T>;
+  if (m < 0 && m != std::numeric_limits<T>::min())
+    m = -m;
+  if (n < 0 && n != std::numeric_limits<T>::min())
+    n = -n;
+  return basic_gcd_(static_cast<Tp>(m), static_cast<Tp>(n));
+}
+
+template <typename Input>
+void do_fuzzy_tests() {
+  std::mt19937 gen(1938);
+  std::uniform_int_distribution<Input> distrib;
+
+  constexpr int nb_rounds = 10000;
+  for (int i = 0; i < nb_rounds; ++i) {
+    Input n = distrib(gen);
+    Input m = distrib(gen);
+    assert(std::gcd(n, m) == basic_gcd(n, m));
+  }
+}
+
+template <typename Input>
+void do_limit_tests() {
+  Input inputs[] = {
+      std::numeric_limits<Input>::min(),
+      std::numeric_limits<Input>::max(),
+      0,
+      1,
+      2,
+      3,
+      4,
+      5,
+      6,
+      7,
+      8,
+      9,
+      10,
+      (Input)-1,
+      (Input)-2,
+      (Input)-3,
+      (Input)-4,
+      (Input)-5,
+      (Input)-6,
+      (Input)-7,
+      (Input)-8,
+      (Input)-9,
+      (Input)-10,
+  };
+
+  for (auto n : inputs) {
+    for (auto m : inputs) {
+      assert(std::gcd(n, m) == basic_gcd(n, m));
+    }
+  }
+}
 
 template <typename Input1, typename Input2 = Input1>
 constexpr bool do_test(int = 0)
@@ -143,5 +206,23 @@ int main(int argc, char**)
     assert(res == 2);
     }
 
-  return 0;
+    do_fuzzy_tests<std::int8_t>();
+    do_fuzzy_tests<std::int16_t>();
+    do_fuzzy_tests<std::int32_t>();
+    do_fuzzy_tests<std::int64_t>();
+    do_fuzzy_tests<std::uint8_t>();
+    do_fuzzy_tests<std::uint16_t>();
+    do_fuzzy_tests<std::uint32_t>();
+    do_fuzzy_tests<std::uint64_t>();
+
+    do_limit_tests<std::int8_t>();
+    do_limit_tests<std::int16_t>();
+    do_limit_tests<std::int32_t>();
+    do_limit_tests<std::int64_t>();
+    do_limit_tests<std::uint8_t>();
+    do_limit_tests<std::uint16_t>();
+    do_limit_tests<std::uint32_t>();
+    do_limit_tests<std::uint64_t>();
+
+    return 0;
 }