Implement QofInt128::gcd and lcm.

src/libqof/qof/qofint128.cpp

@@ -68,7 +68,7 @@ QofInt128::QofInt128 (int64_t upper, int64_t lower, unsigned char flags) :
     m_hi >>= 1;
 }
 
-    QofInt128::QofInt128 (uint64_t upper, uint64_t lower, unsigned char flags) :
+QofInt128::QofInt128 (uint64_t upper, uint64_t lower, unsigned char flags) :
     m_flags {flags}, m_hi {upper},
     m_lo {lower} {}
 
@@ -124,6 +124,55 @@ QofInt128::cmp (const QofInt128& b) const noexcept
     return 0;
 }
 
+/* Knuth 4.5.3 Algo B, recommended by GMP as much faster than Algo A (Euclidean
+ * method).
+ */
+QofInt128
+QofInt128::gcd(QofInt128 b) const noexcept
+{
+    if (b.isZero())
+        return *this;
+    if (isZero())
+        return b;
+
+    if (b.isOverflow() || b.isNan())
+        return b;
+    if (isOverflow() || isNan())
+        return *this;
+
+    QofInt128 a (isNeg() ? -(*this) : *this);
+    if (b.isNeg()) b = -b;
+
+    uint k {};
+    const uint64_t one {1};
+    while (!((a & one) || (b & one))) //B1
+    {
+        a >>= 1;
+        b >>= 1;
+        ++k;
+    }
+    QofInt128 t {a & one ? -b : a}; //B2
+    while (a != b)
+    {
+        while (t && (t & one ^ one)) t >>= 1;  //B3 & B4
+        if (t.isNeg())  //B5
+            b = -t;
+        else
+            a = t;
+        t = a - b;  //B6
+    }
+    return a << k;
+}
+
+/* Since u * v = gcd(u, v) * lcm(u, v), we find lcm by u / gcd * v. */
+
+QofInt128
+QofInt128::lcm(const QofInt128& b) const noexcept
+{
+    auto common = gcd(b);
+    return *this / common * b;
+}
+
 bool
 QofInt128::isNeg () const noexcept
 {

src/libqof/qof/qofint128.hpp

@@ -114,7 +114,7 @@ enum // Values for m_flags
  *
  * @return A QofInt128 having the GCD.
  */
-    QofInt128 gcd (const QofInt128& b) const noexcept;
+    QofInt128 gcd (QofInt128 b) const noexcept;
 /**
  * Computes the Least Common Multiple between the object and parameter
  *

src/libqof/qof/test/gtest-qofint128.cpp

@@ -388,3 +388,26 @@ TEST(qofint128_functions, divide)
     EXPECT_EQ (big, big %= bigger);
     EXPECT_EQ (two, bigger /= big);
 }
+
+TEST(qofint128_functions, GCD)
+{
+    int64_t barg {INT64_C(4878849681579065407)};
+    int64_t sarg {INT64_C(4344522355275710401)};
+    uint64_t uarg {UINT64_C(13567894392130486208)};
+
+    QofInt128 one (INT64_C(1));
+    QofInt128 smallest (sarg);
+    QofInt128 smaller (barg);
+    QofInt128 small (uarg);
+
+    QofInt128 big = smaller * smallest;
+    QofInt128 bigger = small * smaller;
+
+    EXPECT_EQ (smaller, big.gcd(smaller));
+    EXPECT_EQ (smallest, big.gcd(smallest));
+    EXPECT_EQ (small, bigger.gcd(small));
+    EXPECT_EQ (smaller, bigger.gcd(smaller));
+    EXPECT_EQ (one, big.gcd (small));
+    EXPECT_EQ (one, bigger.gcd (smallest));
+    EXPECT_EQ (big, smaller.lcm (smallest));
+}