[libcxx-commits] [libcxx] [libc++] Fix sub-overflow in std::gcd implementation (PR #117984)

[via libcxx-commits]

https://github.com/serge-sans-paille updated https://github.com/llvm/llvm-project/pull/117984

>From 1f04c921953aa5ab36c5cca0464be72b10c29c76 Mon Sep 17 00:00:00 2001
From: serge-sans-paille <sguelton at mozilla.com>
Date: Thu, 28 Nov 2024 10:36:18 +0100
Subject: [PATCH 1/2] [libc++] Fix sub-overflow in std::gcd implementation

Fix #117249
---
 libcxx/include/__numeric/gcd_lcm.h | 22 +++++++++-------------
 1 file changed, 9 insertions(+), 13 deletions(-)

diff --git a/libcxx/include/__numeric/gcd_lcm.h b/libcxx/include/__numeric/gcd_lcm.h
index 9be6cf8516b131..371ca30129a010 100644
--- a/libcxx/include/__numeric/gcd_lcm.h
+++ b/libcxx/include/__numeric/gcd_lcm.h
@@ -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/
   //
   // If power of two divides both numbers, we can push it out.
   // - gcd( 2^x * a, 2^x * b) = 2^x * gcd(a, b)
@@ -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;
 }

>From fcaad75599ba3ee65f4ee4e89bcc7e42a8ad2a1c Mon Sep 17 00:00:00 2001
From: serge-sans-paille <sguelton at mozilla.com>
Date: Thu, 5 Dec 2024 10:36:27 +0100
Subject: [PATCH 2/2] fixup! [libc++] Fix sub-overflow in std::gcd
 implementation

---
 libcxx/include/__numeric/gcd_lcm.h                   |  3 ++-
 .../numeric.ops/numeric.ops.gcd/gcd.pass.cpp         | 12 +-----------
 2 files changed, 3 insertions(+), 12 deletions(-)

diff --git a/libcxx/include/__numeric/gcd_lcm.h b/libcxx/include/__numeric/gcd_lcm.h
index 371ca30129a010..f15f64ea5568d0 100644
--- a/libcxx/include/__numeric/gcd_lcm.h
+++ b/libcxx/include/__numeric/gcd_lcm.h
@@ -55,7 +55,8 @@ template <class _Tp>
 constexpr _LIBCPP_HIDDEN _Tp __gcd(_Tp __a, _Tp __b) {
   static_assert(!is_signed<_Tp>::value, "");
 
-  // From: https://lemire.me/blog/2024/04/13/greatest-common-divisor-the-extended-euclidean-algorithm-and-speed/
+  // Using Binary GCD algorithm https://en.wikipedia.org/wiki/Binary_GCD_algorithm, based on an implementation
+  // from https://lemire.me/blog/2024/04/13/greatest-common-divisor-the-extended-euclidean-algorithm-and-speed/
   //
   // If power of two divides both numbers, we can push it out.
   // - gcd( 2^x * a, 2^x * b) = 2^x * gcd(a, b)
diff --git a/libcxx/test/std/numerics/numeric.ops/numeric.ops.gcd/gcd.pass.cpp b/libcxx/test/std/numerics/numeric.ops/numeric.ops.gcd/gcd.pass.cpp
index 6a9ec1a2ffec24..975b53a763afa4 100644
--- a/libcxx/test/std/numerics/numeric.ops/numeric.ops.gcd/gcd.pass.cpp
+++ b/libcxx/test/std/numerics/numeric.ops/numeric.ops.gcd/gcd.pass.cpp
@@ -27,17 +27,7 @@ constexpr struct {
   int x;
   int y;
   int expect;
-} Cases[] = {
-    {0, 0, 0},
-    {1, 0, 1},
-    {0, 1, 1},
-    {1, 1, 1},
-    {2, 3, 1},
-    {2, 4, 2},
-    {36, 17, 1},
-    {36, 18, 18}
-};
-
+} Cases[] = {{0, 0, 0}, {1, 0, 1}, {0, 1, 1}, {1, 1, 1}, {2, 3, 1}, {2, 4, 2}, {11, 9, 1}, {36, 17, 1}, {36, 18, 18}};
 
 template <typename Input1, typename Input2, typename Output>
 constexpr bool test0(int in1, int in2, int out)