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

[serge-sans-paille]

Summary

This implementation is ~two times faster than current implementation in my setup

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/

[serge-sans-paille]

First time contribution to libc++, I apologize in advance if the style / conformance is not correct.

[serge-sans-paille]

+ constexpr version
+ check for ub in bitwise shift.

[serge-sans-paille]

Still get a

# command output:
*** /home/ssp/sources/llvm-project/libcxx/test/libcxx/transitive_includes/cxx2b.csv
--- /home/ssp/sources/llvm-project/_build/runtimes/runtimes-bins/test/libcxx/Output/transitive_includes.sh.cpp.dir/t.tmp/transitive_includes.csv
***************
*** 432,437 ****
--- 432,438 ----
  numbersversion
  numericcmath
// -*- C++ -*-
  numericcstddef
+ numericinitializer_list
  numericlimits
  numericversion
  optionalcompare

error: command failed with exit status: 1
********************
Failed Tests (1):
  llvm-libc++-shared.cfg.in :: libcxx/transitive_includes.sh.cpp

[cjdb]

Thanks for working on this. I haven't contributed to libc++ in a very long time, so the project's policies may be different to what I'm asking for. Having said that, the following things would be good to do:

• Please add a description to your commit message explaining why this change is valuable. If there's a GitHub issue, please also link to that (just add something like Closes #1234).
• If this is a functional change, please add a test. If this is a non-functional change, please add [NFC] to the commit message.
• It may be worth moving the #ifdef to the body of the function, rather than creating two identical overloads whose interfaces might accidentally diverge over time.

libcxx/include/__numeric/gcd_lcm.h

64

Please fully qualify all the function calls.

[serge-sans-paille]

Take review into account.

[serge-sans-paille]

Fix test

[serge-sans-paille]

Fix ADL

[serge-sans-paille]

@cjdb : should be good now :-)

[serge-sans-paille]

Fix variable naming.

[serge-sans-paille]

rebased + fix intermediate type

[philnik]

Thanks for working on this!
Could you add the benchmark you used to libcxx/benchmarks and give us the exact results you get?

libcxx/include/__numeric/gcd_lcm.h

15

These should be included unconditionally.

67

Same throughout

68

Is there anything C++20-specific other than countr_zero? If no, then I think it would make sense to add a __countr_zero to older language modes and forward the countr_zero implementation.

70

[serge-sans-paille]

Associated benchmark: https://quick-bench.com/q/MpH_u9y6pjBes9p3MFiElANDHXI

[philnik]

Neither quick bench nor I on my local system can reproduce a 4x speedup. In fact, the naive version is significantly faster than your version. Are you sure everything works as expected?

libcxx/include/__bit/countr.h

35 ↗

(On Diff #525484)

This doesn't work pre-C++20.

36 ↗

(On Diff #525484)

libcxx/include/__numeric/gcd_lcm.h

87

Same for other places.

libcxx/test/libcxx/transitive_includes/cxx26.csv

188–208 ↗

(On Diff #525484)

This shouldn't be in here anymore. Did you maybe use an older version to generate this?

[serge-sans-paille]

Updated benchmark that doesn't focus on a single value: https://quick-bench.com/q/XEb2zlF7cpb_gSjUB2f71X4Hbw8

[philnik]

libcxx/include/__bit/countr.h

28 ↗

(On Diff #525994)

Please add the benchmarks under libcxx/benchmarks/.

libcxx/include/__numeric/gcd_lcm.h

59

Could you add a short comment explaining how the algorithm works?

libcxx/include/__numeric/gcd_lcm.h

 // -*- C++ -*-
 //===----------------------------------------------------------------------===//
 //
 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
 // See https://llvm.org/LICENSE.txt for license information.
 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
 //
 //===----------------------------------------------------------------------===//
 
 #ifndef _LIBCPP___NUMERIC_GCD_LCM_H
 #define _LIBCPP___NUMERIC_GCD_LCM_H
 
 #include <__assert>
 #include <__config>
+#if _LIBCPP_STD_VER >= 20
+#  include <__algorithm/min.h>
+#  include <__bit/countr.h>
+#endif
 #include <__type_traits/common_type.h>
 #include <__type_traits/is_integral.h>
 #include <__type_traits/is_same.h>
 #include <__type_traits/is_signed.h>
 #include <__type_traits/make_unsigned.h>
 #include <limits>
 
 #if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER)

 };
 
 template <typename _Result, typename _Source>
 struct __ct_abs<_Result, _Source, false> {
     _LIBCPP_CONSTEXPR _LIBCPP_INLINE_VISIBILITY
     _Result operator()(_Source __t) const noexcept { return __t; }
 };
 
-
 template<class _Tp>
 _LIBCPP_CONSTEXPR _LIBCPP_HIDDEN
 _Tp __gcd(_Tp __m, _Tp __n)
 {
     static_assert((!is_signed<_Tp>::value), "");
+#  if _LIBCPP_STD_VER >= 20
+    if (__m == 0)
+        return __n;
+    if (__n == 0)
+        return __m;
+
+    int __az = _VSTD::countr_zero(__m);
+    int __bz = _VSTD::countr_zero(__n);
+    int __shift = _VSTD::min(__az, __bz);
+    __m >>= __az, __n >>= __bz;
+    while (__m != 0) {
+        int diff = __m - __n;
+        __n = _VSTD::min(__m, __n);
+        __m = diff >= 0 ? diff : -diff;
+        int __m_shift = _VSTD::min(countr_zero(__m), _VSTD::numeric_limits<_Tp>::digits - 1);
+        __m >>= __m_shift;
+    }
+    return __n << __shift;
+#  else
     return __n == 0 ? __m : _VSTD::__gcd<_Tp>(__n, __m % __n);
+#  endif
 }
 
 template<class _Tp, class _Up>
 _LIBCPP_CONSTEXPR _LIBCPP_INLINE_VISIBILITY
 common_type_t<_Tp,_Up>
 gcd(_Tp __m, _Up __n)
 {
     static_assert((is_integral<_Tp>::value && is_integral<_Up>::value), "Arguments to gcd must be integer types");
     static_assert((!is_same<__remove_cv_t<_Tp>, bool>::value), "First argument to gcd cannot be bool" );
     static_assert((!is_same<__remove_cv_t<_Up>, bool>::value), "Second argument to gcd cannot be bool" );
     using _Rp = common_type_t<_Tp,_Up>;
     using _Wp = make_unsigned_t<_Rp>;
     return static_cast<_Rp>(_VSTD::__gcd(
         static_cast<_Wp>(__ct_abs<_Rp, _Tp>()(__m)),