Skip to content
Merged
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 9 additions & 13 deletions libcxx/include/__numeric/gcd_lcm.h
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,7 @@ template <class _Tp>
constexpr _LIBCPP_HIDDEN _Tp __gcd(_Tp __a, _Tp __b) {
static_assert(!is_signed<_Tp>::value, "");

// From: https://lemire.me/blog/2013/12/26/fastest-way-to-compute-the-greatest-common-divisor
// From: https://lemire.me/blog/2024/04/13/greatest-common-divisor-the-extended-euclidean-algorithm-and-speed/
Copy link
Collaborator

@hiraditya hiraditya Dec 2, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I believe this is Stein's algorithm (Can't find the link currently but this has a reference on https://aha.stanford.edu/sites/g/files/sbiybj20066/files/media/file/aha_050422_sreedhar_crypto_xgcd_1.pdf#page=14). Maybe we should write that instead of linking to a blogpost? But I'm not sure what the policy is about adding a link in the code. Can we add the blogpost link in a commit message instead? There is a high chance blogpost may be dead in a few years.

//
// If power of two divides both numbers, we can push it out.
// - gcd( 2^x * a, 2^x * b) = 2^x * gcd(a, b)
Expand All @@ -76,21 +76,17 @@ constexpr _LIBCPP_HIDDEN _Tp __gcd(_Tp __a, _Tp __b) {
if (__a == 0)
return __b;

int __az = std::__countr_zero(__a);
int __bz = std::__countr_zero(__b);
int __shift = std::min(__az, __bz);
__a >>= __az;
__b >>= __bz;
_Tp __c = __a | __b;
int __shift = std::__countr_zero(__c);
__a >>= std::__countr_zero(__a);
do {
_Tp __diff = __a - __b;
if (__a > __b) {
__a = __b;
__b = __diff;
_Tp __t = __b >> std::__countr_zero(__b);
if (__a > __t) {
__b = __a - __t;
__a = __t;
} else {
__b = __b - __a;
__b = __t - __a;
}
if (__diff != 0)
__b >>= std::__countr_zero(__diff);
} while (__b != 0);
return __a << __shift;
}
Expand Down