More carry calculation fixing (#835)

* Add carry calculation tests

* Fix carry calculation

* Minor simplification

* Test both ways, now that we are here

* Add helpful remark

* Avoid a single multiplication instruction per product (=one for `a*b`, and another for `c*d`)

CPP/Clipper2Lib/include/clipper2/clipper.core.h

@@ -664,42 +664,51 @@ namespace Clipper2Lib
     return (x > 0) - (x < 0); 
   }
 
-  inline uint64_t CalcOverflowCarry(uint64_t a, uint64_t b) // #834
+  struct MultiplicationResult
   {
-    // given aLo = (a & 0xFFFFFFFF) and
-    // aHi = (a & 0xFFFFFFFF00000000) and similarly with b, then
-    // a * b == (aHi + aLo) * (bHi + bLo)
-    // a * b == (aHi * bHi) + (aHi * bLo) + (aLo * bHi) + (aLo * bLo)
-
-    const uint64_t aLo = a & 0xFFFFFFFF;
-    const uint64_t aHi = a >> 32; // this avoids repeating shifts
-    const uint64_t bLo = b & 0xFFFFFFFF;
-    const uint64_t bHi = b >> 32; 
-    // integer overflow of multiplying the unsigned 64bits a and b ==>
-    return aHi * bHi + ((aHi * bLo) >> 32) + ((bHi * aLo) >> 32);
+    const uint64_t result = 0;
+    const uint64_t carry = 0;
+
+    bool operator==(const MultiplicationResult& other) const
+    {
+      return result == other.result && carry == other.carry;
+    };
+  };
+
+  inline MultiplicationResult Multiply(uint64_t a, uint64_t b) // #834
+  {
+    const auto lo = [](uint64_t x) { return x & 0xFFFFFFFF; };
+    const auto hi = [](uint64_t x) { return x >> 32; };
+
+    // https://stackoverflow.com/a/1815371/1158913
+    const uint64_t x1 = lo(a) * lo(b);
+    const uint64_t x2 = hi(a) * lo(b) + hi(x1);
+    const uint64_t x3 = lo(a) * hi(b) + lo(x2);
+    const uint64_t x4 = hi(a) * hi(b) + hi(x2) + hi(x3);
+
+    const uint64_t result = lo(x3) << 32 | lo(x1);
+    const uint64_t carry  = x4;
+
+    return { result, carry };
   }
 
   // returns true if (and only if) a * b == c * d
   inline bool ProductsAreEqual(int64_t a, int64_t b, int64_t c, int64_t d)
   {
-    // nb: unsigned values will be needed for CalcOverflowCarry()
+    // nb: unsigned values needed for calculating overflow carry
     const auto abs_a = static_cast<uint64_t>(std::abs(a));
     const auto abs_b = static_cast<uint64_t>(std::abs(b));
     const auto abs_c = static_cast<uint64_t>(std::abs(c));
     const auto abs_d = static_cast<uint64_t>(std::abs(d));
 
-    // the multiplications here can potentially overflow, but
-    // any overflows will be compared using CalcOverflowCarry()
-    const auto abs_ab = abs_a * abs_b;
-    const auto abs_cd = abs_c * abs_d;
+    const auto abs_ab = Multiply(abs_a, abs_b);
+    const auto abs_cd = Multiply(abs_c, abs_d);
 
     // nb: it's important to differentiate 0 values here from other values
     const auto sign_ab = TriSign(a) * TriSign(b);
     const auto sign_cd = TriSign(c) * TriSign(d);
 
-    const auto carry_ab = CalcOverflowCarry(abs_a, abs_b);
-    const auto carry_cd = CalcOverflowCarry(abs_c, abs_d);
-    return abs_ab == abs_cd && sign_ab == sign_cd && carry_ab == carry_cd;
+    return abs_ab == abs_cd && sign_ab == sign_cd;
   }
 
   template <typename T>

CPP/Tests/TestIsCollinear.cpp

@@ -1,6 +1,29 @@
 #include <gtest/gtest.h>
 #include "clipper2/clipper.h"
 
+TEST(Clipper2Tests, TestCarryCalculation) {
+  EXPECT_EQ(Clipper2Lib::Multiply(0x51eaed81157de061, 0x3a271fb2745b6fe9).carry, 0x129bbebdfae0464e);
+  EXPECT_EQ(Clipper2Lib::Multiply(0x3a271fb2745b6fe9, 0x51eaed81157de061).carry, 0x129bbebdfae0464e);
+  EXPECT_EQ(Clipper2Lib::Multiply(0xc2055706a62883fa, 0x26c78bc79c2322cc).carry, 0x1d640701d192519b);
+  EXPECT_EQ(Clipper2Lib::Multiply(0x26c78bc79c2322cc, 0xc2055706a62883fa).carry, 0x1d640701d192519b);
+  EXPECT_EQ(Clipper2Lib::Multiply(0x874ddae32094b0de, 0x9b1559a06fdf83e0).carry, 0x51f76c49563e5bfe);
+  EXPECT_EQ(Clipper2Lib::Multiply(0x9b1559a06fdf83e0, 0x874ddae32094b0de).carry, 0x51f76c49563e5bfe);
+  EXPECT_EQ(Clipper2Lib::Multiply(0x81fb3ad3636ca900, 0x239c000a982a8da4).carry, 0x12148e28207b83a3);
+  EXPECT_EQ(Clipper2Lib::Multiply(0x239c000a982a8da4, 0x81fb3ad3636ca900).carry, 0x12148e28207b83a3);
+  EXPECT_EQ(Clipper2Lib::Multiply(0x4be0b4c5d2725c44, 0x990cd6db34a04c30).carry, 0x2d5d1a4183fd6165);
+  EXPECT_EQ(Clipper2Lib::Multiply(0x990cd6db34a04c30, 0x4be0b4c5d2725c44).carry, 0x2d5d1a4183fd6165);
+  EXPECT_EQ(Clipper2Lib::Multiply(0x978ec0c0433c01f6, 0x2df03d097966b536).carry, 0x1b3251d91fe272a5);
+  EXPECT_EQ(Clipper2Lib::Multiply(0x2df03d097966b536, 0x978ec0c0433c01f6).carry, 0x1b3251d91fe272a5);
+  EXPECT_EQ(Clipper2Lib::Multiply(0x49c5cbbcfd716344, 0xc489e3b34b007ad3).carry, 0x38a32c74c8c191a4);
+  EXPECT_EQ(Clipper2Lib::Multiply(0xc489e3b34b007ad3, 0x49c5cbbcfd716344).carry, 0x38a32c74c8c191a4);
+  EXPECT_EQ(Clipper2Lib::Multiply(0xd3361cdbeed655d5, 0x1240da41e324953a).carry, 0x0f0f4fa11e7e8f2a);
+  EXPECT_EQ(Clipper2Lib::Multiply(0x1240da41e324953a, 0xd3361cdbeed655d5).carry, 0x0f0f4fa11e7e8f2a);
+  EXPECT_EQ(Clipper2Lib::Multiply(0x51b854f8e71b0ae0, 0x6f8d438aae530af5).carry, 0x239c04ee3c8cc248);
+  EXPECT_EQ(Clipper2Lib::Multiply(0x6f8d438aae530af5, 0x51b854f8e71b0ae0).carry, 0x239c04ee3c8cc248);
+  EXPECT_EQ(Clipper2Lib::Multiply(0xbbecf7dbc6147480, 0xbb0f73d0f82e2236).carry, 0x895170f4e9a216a7);
+  EXPECT_EQ(Clipper2Lib::Multiply(0xbb0f73d0f82e2236, 0xbbecf7dbc6147480).carry, 0x895170f4e9a216a7);
+}
+
 TEST(Clipper2Tests, TestIsCollinear) {
   // a large integer not representable by double
   const int64_t i = 9007199254740993;