131 lines · c
1// -*- C++ -*-2//===----------------------------------------------------------------------===//3//4// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.5// See https://llvm.org/LICENSE.txt for license information.6// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception7//8//===----------------------------------------------------------------------===//9 10#ifndef _LIBCPP___NUMERIC_GCD_LCM_H11#define _LIBCPP___NUMERIC_GCD_LCM_H12 13#include <__assert>14#include <__bit/countr.h>15#include <__config>16#include <__memory/addressof.h>17#include <__type_traits/common_type.h>18#include <__type_traits/is_integral.h>19#include <__type_traits/is_same.h>20#include <__type_traits/is_signed.h>21#include <__type_traits/make_unsigned.h>22#include <__type_traits/remove_cv.h>23#include <limits>24 25#if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER)26# pragma GCC system_header27#endif28 29_LIBCPP_PUSH_MACROS30#include <__undef_macros>31 32_LIBCPP_BEGIN_NAMESPACE_STD33 34#if _LIBCPP_STD_VER >= 1735 36template <typename _Result, typename _Source, bool _IsSigned = is_signed<_Source>::value>37struct __ct_abs;38 39template <typename _Result, typename _Source>40struct __ct_abs<_Result, _Source, true> {41 constexpr _LIBCPP_HIDE_FROM_ABI _Result operator()(_Source __t) const noexcept {42 if (__t >= 0)43 return __t;44 if (__t == numeric_limits<_Source>::min())45 return -static_cast<_Result>(__t);46 return -__t;47 }48};49 50template <typename _Result, typename _Source>51struct __ct_abs<_Result, _Source, false> {52 constexpr _LIBCPP_HIDE_FROM_ABI _Result operator()(_Source __t) const noexcept { return __t; }53};54 55template <class _Tp>56constexpr _LIBCPP_HIDDEN _Tp __gcd(_Tp __a, _Tp __b) {57 static_assert(!is_signed<_Tp>::value, "");58 59 // Using Binary GCD algorithm https://en.wikipedia.org/wiki/Binary_GCD_algorithm, based on an implementation60 // from https://lemire.me/blog/2024/04/13/greatest-common-divisor-the-extended-euclidean-algorithm-and-speed/61 //62 // If power of two divides both numbers, we can push it out.63 // - gcd( 2^x * a, 2^x * b) = 2^x * gcd(a, b)64 //65 // If and only if exactly one number is even, we can divide that number by that power.66 // - if a, b are odd, then gcd(2^x * a, b) = gcd(a, b)67 //68 // And standard gcd algorithm where instead of modulo, minus is used.69 70 if (__a < __b) {71 _Tp __tmp = __b;72 __b = __a;73 __a = __tmp;74 }75 if (__b == 0)76 return __a;77 __a %= __b; // Make both argument of the same size, and early result in the easy case.78 if (__a == 0)79 return __b;80 81 _Tp __c = __a | __b;82 int __shift = std::__countr_zero(__c);83 __a >>= std::__countr_zero(__a);84 do {85 _Tp __t = __b >> std::__countr_zero(__b);86 if (__a > __t) {87 __b = __a - __t;88 __a = __t;89 } else {90 __b = __t - __a;91 }92 } while (__b != 0);93 return __a << __shift;94}95 96template <class _Tp, class _Up>97constexpr _LIBCPP_HIDE_FROM_ABI common_type_t<_Tp, _Up> gcd(_Tp __m, _Up __n) {98 static_assert(is_integral<_Tp>::value && is_integral<_Up>::value, "Arguments to gcd must be integer types");99 static_assert(!is_same<__remove_cv_t<_Tp>, bool>::value, "First argument to gcd cannot be bool");100 static_assert(!is_same<__remove_cv_t<_Up>, bool>::value, "Second argument to gcd cannot be bool");101 using _Rp = common_type_t<_Tp, _Up>;102 using _Wp = make_unsigned_t<_Rp>;103 return static_cast<_Rp>(104 std::__gcd(static_cast<_Wp>(__ct_abs<_Rp, _Tp>()(__m)), static_cast<_Wp>(__ct_abs<_Rp, _Up>()(__n))));105}106 107template <class _Tp, class _Up>108constexpr _LIBCPP_HIDE_FROM_ABI common_type_t<_Tp, _Up> lcm(_Tp __m, _Up __n) {109 static_assert(is_integral<_Tp>::value && is_integral<_Up>::value, "Arguments to lcm must be integer types");110 static_assert(!is_same<__remove_cv_t<_Tp>, bool>::value, "First argument to lcm cannot be bool");111 static_assert(!is_same<__remove_cv_t<_Up>, bool>::value, "Second argument to lcm cannot be bool");112 if (__m == 0 || __n == 0)113 return 0;114 115 using _Rp = common_type_t<_Tp, _Up>;116 _Rp __val1 = __ct_abs<_Rp, _Tp>()(__m) / std::gcd(__m, __n);117 _Rp __val2 = __ct_abs<_Rp, _Up>()(__n);118 _Rp __res;119 [[maybe_unused]] bool __overflow = __builtin_mul_overflow(__val1, __val2, std::addressof(__res));120 _LIBCPP_ASSERT_ARGUMENT_WITHIN_DOMAIN(!__overflow, "Overflow in lcm");121 return __res;122}123 124#endif // _LIBCPP_STD_VER >= 17125 126_LIBCPP_END_NAMESPACE_STD127 128_LIBCPP_POP_MACROS129 130#endif // _LIBCPP___NUMERIC_GCD_LCM_H131