v1/__numeric/gcd_lcm.h

0001 // -*- C++ -*-
0002 //===----------------------------------------------------------------------===//
0003 //
0004 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
0005 // See https://llvm.org/LICENSE.txt for license information.
0006 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
0007 //
0008 //===----------------------------------------------------------------------===//
0009
0010 #ifndef _LIBCPP___NUMERIC_GCD_LCM_H
0011 #define _LIBCPP___NUMERIC_GCD_LCM_H
0012
0013 #include <__algorithm/min.h>
0014 #include <__assert>
0015 #include <__bit/countr.h>
0016 #include <__config>
0017 #include <__type_traits/common_type.h>
0018 #include <__type_traits/is_integral.h>
0019 #include <__type_traits/is_same.h>
0020 #include <__type_traits/is_signed.h>
0021 #include <__type_traits/make_unsigned.h>
0022 #include <limits>
0023
0024 #if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER)
0025 #  pragma GCC system_header
0026 #endif
0027
0028 _LIBCPP_PUSH_MACROS
0029 #include <__undef_macros>
0030
0031 _LIBCPP_BEGIN_NAMESPACE_STD
0032
0033 #if _LIBCPP_STD_VER >= 17
0034
0035 template <typename _Result, typename _Source, bool _IsSigned = is_signed<_Source>::value>
0036 struct __ct_abs;
0037
0038 template <typename _Result, typename _Source>
0039 struct __ct_abs<_Result, _Source, true> {
0040   constexpr _LIBCPP_HIDE_FROM_ABI _Result operator()(_Source __t) const noexcept {
0041     if (__t >= 0)
0042       return __t;
0043     if (__t == numeric_limits<_Source>::min())
0044       return -static_cast<_Result>(__t);
0045     return -__t;
0046   }
0047 };
0048
0049 template <typename _Result, typename _Source>
0050 struct __ct_abs<_Result, _Source, false> {
0051   constexpr _LIBCPP_HIDE_FROM_ABI _Result operator()(_Source __t) const noexcept { return __t; }
0052 };
0053
0054 template <class _Tp>
0055 constexpr _LIBCPP_HIDDEN _Tp __gcd(_Tp __a, _Tp __b) {
0056   static_assert(!is_signed<_Tp>::value, "");
0057
0058   // Using Binary GCD algorithm https://en.wikipedia.org/wiki/Binary_GCD_algorithm, based on an implementation
0059   // from https://lemire.me/blog/2024/04/13/greatest-common-divisor-the-extended-euclidean-algorithm-and-speed/
0060   //
0061   // If power of two divides both numbers, we can push it out.
0062   // - gcd( 2^x * a, 2^x * b) = 2^x * gcd(a, b)
0063   //
0064   // If and only if exactly one number is even, we can divide that number by that power.
0065   // - if a, b are odd, then gcd(2^x * a, b) = gcd(a, b)
0066   //
0067   // And standard gcd algorithm where instead of modulo, minus is used.
0068
0069   if (__a < __b) {
0070     _Tp __tmp = __b;
0071     __b       = __a;
0072     __a       = __tmp;
0073   }
0074   if (__b == 0)
0075     return __a;
0076   __a %= __b; // Make both argument of the same size, and early result in the easy case.
0077   if (__a == 0)
0078     return __b;
0079
0080   _Tp __c     = __a | __b;
0081   int __shift = std::__countr_zero(__c);
0082   __a >>= std::__countr_zero(__a);
0083   do {
0084     _Tp __t = __b >> std::__countr_zero(__b);
0085     if (__a > __t) {
0086       __b = __a - __t;
0087       __a = __t;
0088     } else {
0089       __b = __t - __a;
0090     }
0091   } while (__b != 0);
0092   return __a << __shift;
0093 }
0094
0095 template <class _Tp, class _Up>
0096 constexpr _LIBCPP_HIDE_FROM_ABI common_type_t<_Tp, _Up> gcd(_Tp __m, _Up __n) {
0097   static_assert(is_integral<_Tp>::value && is_integral<_Up>::value, "Arguments to gcd must be integer types");
0098   static_assert(!is_same<__remove_cv_t<_Tp>, bool>::value, "First argument to gcd cannot be bool");
0099   static_assert(!is_same<__remove_cv_t<_Up>, bool>::value, "Second argument to gcd cannot be bool");
0100   using _Rp = common_type_t<_Tp, _Up>;
0101   using _Wp = make_unsigned_t<_Rp>;
0102   return static_cast<_Rp>(
0103       std::__gcd(static_cast<_Wp>(__ct_abs<_Rp, _Tp>()(__m)), static_cast<_Wp>(__ct_abs<_Rp, _Up>()(__n))));
0104 }
0105
0106 template <class _Tp, class _Up>
0107 constexpr _LIBCPP_HIDE_FROM_ABI common_type_t<_Tp, _Up> lcm(_Tp __m, _Up __n) {
0108   static_assert(is_integral<_Tp>::value && is_integral<_Up>::value, "Arguments to lcm must be integer types");
0109   static_assert(!is_same<__remove_cv_t<_Tp>, bool>::value, "First argument to lcm cannot be bool");
0110   static_assert(!is_same<__remove_cv_t<_Up>, bool>::value, "Second argument to lcm cannot be bool");
0111   if (__m == 0 || __n == 0)
0112     return 0;
0113
0114   using _Rp  = common_type_t<_Tp, _Up>;
0115   _Rp __val1 = __ct_abs<_Rp, _Tp>()(__m) / std::gcd(__m, __n);
0116   _Rp __val2 = __ct_abs<_Rp, _Up>()(__n);
0117   _Rp __res;
0118   [[maybe_unused]] bool __overflow = __builtin_mul_overflow(__val1, __val2, &__res);
0119   _LIBCPP_ASSERT_ARGUMENT_WITHIN_DOMAIN(!__overflow, "Overflow in lcm");
0120   return __res;
0121 }
0122
0123 #endif // _LIBCPP_STD_VER >= 17
0124
0125 _LIBCPP_END_NAMESPACE_STD
0126
0127 _LIBCPP_POP_MACROS
0128
0129 #endif // _LIBCPP___NUMERIC_GCD_LCM_H