[libc][math] Implement full multiplication and quick_mul_hi for UInt class.

[lntue]

Summary

Implement full multiplication UInt<A> * UInt<B> -> UInt<A + B> and
quick_mul_hi that returns the higher half of the product UInt<A> * UInt<A>.
These 2 functions will be used for dyadic floating point class.

[sivachandra]

libc/src/__support/UInt.h

219

Can you add more commentary here explaining how this is different from operator*? IIUC, operator* returns the LS Bits of the result, while this is will return the MS Bits of the result with an additional condition that the LS-Bit of the result need not be the same as that in the MS Bits of the result produced by ful_mul.

libc/test/src/__support/uint_test.cpp

510

Please fix.

libc/utils/UnitTest/LibcTest.cpp

294

While OK for this change, I think this is going beyond manageable limits. We should refactor the test infrastructure with some urgency so that we can put these specializations next to the special types.

[lntue]

Add comments about error bounds and efficiency of quick_mul_hi + More tests.

[sivachandra]

libc/src/__support/UInt.h

219

Add an example of where this kind of multiplication will be useful.

libc/src/__support/UInt.h

if constexpr (WordCount == 1) {
         partial_sum.val[0] = partial_sum.val[1];
         partial_sum.val[1] = carry;
         carry = 0;
       }
       return result;
     }
   }
 
+  // Return the full product.
+  template <size_t OtherBits>
+  constexpr UInt<Bits + OtherBits> ful_mul(const UInt<OtherBits> &other) const {
+    UInt<Bits + OtherBits> result(0);
+    UInt<128> partial_sum(0);
+    uint64_t carry = 0;
+    constexpr size_t OtherWordCount = UInt<OtherBits>::WordCount;
+    for (size_t i = 0; i <= WordCount + OtherWordCount - 2; ++i) {
+      const size_t lower_idx = i < OtherWordCount ? 0 : i - OtherWordCount + 1;
+      const size_t upper_idx = i < WordCount ? i : WordCount - 1;
+      for (size_t j = lower_idx; j <= upper_idx; ++j) {
+        NumberPair<uint64_t> prod = full_mul(val[j], other.val[i - j]);
+        UInt<128> tmp({prod.lo, prod.hi});
+        carry += partial_sum.add(tmp);
+      }
+      result.val[i] = partial_sum.val[0];
+      partial_sum.val[0] = partial_sum.val[1];
+      partial_sum.val[1] = carry;
+      carry = 0;
+    }
+    result.val[WordCount + OtherWordCount - 1] = partial_sum.val[0];
+    return result;
+  }
+
+  // Fast hi part of the full product.  The normal product `operator*` returns
+  // `Bits` least significant bits of the full product, while this function will
+  // approximate `Bits` most significant bits of the full product with errors
+  // bounded by:
+  //   0 <= (a.full_mul(b) >> Bits) - a.quick_mul_hi(b)) <= WordCount - 1.
+  //
+  // An example usage of this is to quickly (but less accurately) compute the
+  // product of (normalized) mantissas of floating point numbers:
+  //   (mant_1, mant_2) -> quick_mul_hi -> normalize leading bit
+  // is much more efficient than:
+  //   (mant_1, mant_2) -> ful_mul -> normalize leading bit
+  //                    -> convert back to same Bits width by shifting/rounding,
+  // especially for higher precisions.
+  //
+  // Performance summary:
+  //   Number of 64-bit x 64-bit -> 128-bit multiplications performed.
+  //   Bits  WordCount  ful_mul  quick_mul_hi  Error bound
+  //    128      2         4           3            1
+  //    196      3         9           6            2
+  //    256      4        16          10            3
+  //    512      8        64          36            7
+  constexpr UInt<Bits> quick_mul_hi(const UInt<Bits> &other) const {
+    UInt<Bits> result(0);
+    UInt<128> partial_sum(0);
+    uint64_t carry = 0;
+    // First round of accumulation for those at WordCount - 1 in the full
+    // product.
+    for (size_t i = 0; i < WordCount; ++i) {
+      NumberPair<uint64_t> prod =
+          full_mul(val[i], other.val[WordCount - 1 - i]);
+      UInt<128> tmp({prod.lo, prod.hi});
+      carry += partial_sum.add(tmp);
+    }
+    for (size_t i = WordCount; i < 2 * WordCount - 1; ++i) {
+      partial_sum.val[0] = partial_sum.val[1];
+      partial_sum.val[1] = carry;
+      carry = 0;
+      for (size_t j = i - WordCount + 1; j < WordCount; ++j) {
+        NumberPair<uint64_t> prod = full_mul(val[j], other.val[i - j]);
+        UInt<128> tmp({prod.lo, prod.hi});
+        carry += partial_sum.add(tmp);
+      }
+      result.val[i - WordCount] = partial_sum.val[0];
+    }
+    result.val[WordCount - 1] = partial_sum.val[1];
+    return result;
+  }
+
   // pow takes a power and sets this to its starting value to that power. Zero
   // to the zeroth power returns 1.
   constexpr void pow_n(uint64_t power) {
     UInt<Bits> result = 1;
     UInt<Bits> cur_power = *this;
 
     while (power > 0) {
       if ((power % 2) > 0) {

libc/test/src/__support/CMakeLists.txt

SUITE
     libc_support_unittests
   SRCS
     arg_list_test.cpp
   DEPENDS
     libc.src.__support.arg_list
 )
 
 add_libc_unittest(
-  uint128_test
+  uint_test
   SUITE
     libc_support_unittests
   SRCS
-    uint128_test.cpp
+    uint_test.cpp
   DEPENDS
     libc.src.__support.uint
     libc.src.__support.CPP.optional
 )
 
 add_libc_unittest(
   fixedvector_test
   SUITE

libc/test/src/__support/uint_test.cpp

-//===-- Unittests for the 128 bit integer class ---------------------------===//
+//===-- Unittests for the UInt integer class ------------------------------===//
 //
 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
 // See https://llvm.org/LICENSE.txt for license information.
 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
 //
 //===----------------------------------------------------------------------===//
 
 #include "src/__support/CPP/optional.h"
 #include "src/__support/UInt.h"
 
 #include "utils/UnitTest/Test.h"
 
 // We want to test __llvm_libc::cpp::UInt<128> explicitly. So, for convenience,
 // we use a sugar which does not conflict with the UInt128 type which can
 // resolve to __uint128_t if the platform has it.
 using LL_UInt128 = __llvm_libc::cpp::UInt<128>;
 using LL_UInt192 = __llvm_libc::cpp::UInt<192>;
 using LL_UInt256 = __llvm_libc::cpp::UInt<256>;
+using LL_UInt320 = __llvm_libc::cpp::UInt<320>;
+using LL_UInt512 = __llvm_libc::cpp::UInt<512>;
+using LL_UInt1024 = __llvm_libc::cpp::UInt<1024>;
 
-TEST(LlvmLibcUInt128ClassTest, BasicInit) {
+TEST(LlvmLibcUIntClassTest, BasicInit) {
   LL_UInt128 empty;
   LL_UInt128 half_val(12345);
   LL_UInt128 full_val({12345, 67890});
   ASSERT_TRUE(half_val != full_val);
 }
 
-TEST(LlvmLibcUInt128ClassTest, AdditionTests) {
+TEST(LlvmLibcUIntClassTest, AdditionTests) {
   LL_UInt128 val1(12345);
   LL_UInt128 val2(54321);
   LL_UInt128 result1(66666);
   EXPECT_EQ(val1 + val2, result1);
   EXPECT_EQ((val1 + val2), (val2 + val1)); // addition is commutative
 
   // Test overflow
   LL_UInt128 val3({0xf000000000000001, 0});
TEST(LlvmLibcUIntClassTest, AdditionTests) {
   LL_UInt256 val10({0x1111222233334444, 0xaaaabbbbccccdddd, 0x1111222233334444,
                     0xaaaabbbbccccdddd});
   LL_UInt256 result5({0x302f3f3e4e4d5d5c, 0x9c9dafb0c2c3d5d5,
                       0x12346789bcdf1234, 0xa987765443210fed});
   EXPECT_EQ(val9 + val10, result5);
   EXPECT_EQ(val9 + val10, val10 + val9);
 }
 
-TEST(LlvmLibcUInt128ClassTest, SubtractionTests) {
+TEST(LlvmLibcUIntClassTest, SubtractionTests) {
   LL_UInt128 val1(12345);
   LL_UInt128 val2(54321);
   LL_UInt128 result1({0xffffffffffff5c08, 0xffffffffffffffff});
   LL_UInt128 result2(0xa3f8);
   EXPECT_EQ(val1 - val2, result1);
   EXPECT_EQ(val1, val2 + result1);
   EXPECT_EQ(val2 - val1, result2);
   EXPECT_EQ(val2, val1 + result2);
TEST(LlvmLibcUIntClassTest, SubtractionTests) {
   LL_UInt128 result5({0xf0122345567889ab, 0x5431fedca9875432});
   LL_UInt128 result6({0x0feddcbaa9877655, 0xabce01235678abcd});
   EXPECT_EQ(val5 - val6, result5);
   EXPECT_EQ(val5, val6 + result5);
   EXPECT_EQ(val6 - val5, result6);
   EXPECT_EQ(val6, val5 + result6);
 }
 
-TEST(LlvmLibcUInt128ClassTest, MultiplicationTests) {
+TEST(LlvmLibcUIntClassTest, MultiplicationTests) {
   LL_UInt128 val1({5, 0});
   LL_UInt128 val2({10, 0});
   LL_UInt128 result1({50, 0});
   EXPECT_EQ((val1 * val2), result1);
   EXPECT_EQ((val1 * val2), (val2 * val1)); // multiplication is commutative
 
   // Check that the multiplication works accross the whole number
   LL_UInt128 val3({0xf, 0});
TEST(LlvmLibcUIntClassTest, MultiplicationTests) {
   LL_UInt256 val14({0xffffffffffffffff, 0xffffffffffffffff, 0x3792F412CB06794D,
                     0xCDB02555653131B6});
   LL_UInt256 result7({0x0000000000000001, 0xfe289dbdd36b9a6f,
                       0x291de4c71d5f646c, 0xfd37221cb06d4978});
   EXPECT_EQ((val13 * val14), result7);
   EXPECT_EQ((val13 * val14), (val14 * val13));
 }
 
-TEST(LlvmLibcUInt128ClassTest, DivisionTests) {
+TEST(LlvmLibcUIntClassTest, DivisionTests) {
   LL_UInt128 val1({10, 0});
   LL_UInt128 val2({5, 0});
   LL_UInt128 result1({2, 0});
   EXPECT_EQ((val1 / val2), result1);
   EXPECT_EQ((val1 / result1), val2);
 
   // Check that the division works accross the whole number
   LL_UInt128 val3({0xffffffffffffffff, 0xffffffffffffffff});
TEST(LlvmLibcUIntClassTest, DivisionTests) {
   EXPECT_EQ((val11 / val12), result6);
 
   // Make sure that division handles dividing by zero correctly.
   LL_UInt128 val13({1234, 0});
   LL_UInt128 val14({0, 0});
   EXPECT_FALSE(val13.div(val14).has_value());
 }
 
-TEST(LlvmLibcUInt128ClassTest, ModuloTests) {
+TEST(LlvmLibcUIntClassTest, ModuloTests) {
   LL_UInt128 val1({10, 0});
   LL_UInt128 val2({5, 0});
   LL_UInt128 result1({0, 0});
   EXPECT_EQ((val1 % val2), result1);
 
   LL_UInt128 val3({101, 0});
   LL_UInt128 val4({10, 0});
   LL_UInt128 result2({1, 0});
TEST(LlvmLibcUIntClassTest, ModuloTests) {
   EXPECT_EQ((val15 % val16), result8);
 
   LL_UInt128 val17({5076944270305263619, 54210108624}); // (10 ^ 30) + 3
   LL_UInt128 val18({10, 0});
   LL_UInt128 result9({3, 0});
   EXPECT_EQ((val17 % val18), result9);
 }
 
-TEST(LlvmLibcUInt128ClassTest, PowerTests) {
+TEST(LlvmLibcUIntClassTest, PowerTests) {
   LL_UInt128 val1({10, 0});
   val1.pow_n(30);
   LL_UInt128 result1({5076944270305263616, 54210108624}); // (10 ^ 30)
   EXPECT_EQ(val1, result1);
 
   LL_UInt128 val2({1, 0});
   val2.pow_n(10);
   LL_UInt128 result2({1, 0});
for (size_t i = 0; i < 128; ++i) {
     val_pow_two = 2;
     val_pow_two.pow_n(i);
     result_pow_two = 1;
     result_pow_two = result_pow_two << i;
     EXPECT_EQ(val_pow_two, result_pow_two);
   }
 }
 
-TEST(LlvmLibcUInt128ClassTest, ShiftLeftTests) {
+TEST(LlvmLibcUIntClassTest, ShiftLeftTests) {
   LL_UInt128 val1(0x0123456789abcdef);
   LL_UInt128 result1(0x123456789abcdef0);
   EXPECT_EQ((val1 << 4), result1);
 
   LL_UInt128 val2({0x13579bdf02468ace, 0x123456789abcdef0});
   LL_UInt128 result2({0x02468ace00000000, 0x9abcdef013579bdf});
   EXPECT_EQ((val2 << 32), result2);
   LL_UInt128 val22 = val2;
TEST(LlvmLibcUIntClassTest, ShiftLeftTests) {
   LL_UInt128 result5({0, 0x2468ace000000000});
   EXPECT_EQ((val2 << 100), result5);
 
   LL_UInt128 result6({0, 0});
   EXPECT_EQ((val2 << 128), result6);
   EXPECT_EQ((val2 << 256), result6);
 }
 
-TEST(LlvmLibcUInt128ClassTest, ShiftRightTests) {
+TEST(LlvmLibcUIntClassTest, ShiftRightTests) {
   LL_UInt128 val1(0x0123456789abcdef);
   LL_UInt128 result1(0x00123456789abcde);
   EXPECT_EQ((val1 >> 4), result1);
 
   LL_UInt128 val2({0x13579bdf02468ace, 0x123456789abcdef0});
   LL_UInt128 result2({0x9abcdef013579bdf, 0x0000000012345678});
   EXPECT_EQ((val2 >> 32), result2);
   LL_UInt128 val22 = val2;
TEST(LlvmLibcUIntClassTest, ShiftRightTests) {
   LL_UInt128 result5({0x0000000001234567, 0});
   EXPECT_EQ((val2 >> 100), result5);
 
   LL_UInt128 result6({0, 0});
   EXPECT_EQ((val2 >> 128), result6);
   EXPECT_EQ((val2 >> 256), result6);
 }
 
-TEST(LlvmLibcUInt128ClassTest, AndTests) {
+TEST(LlvmLibcUIntClassTest, AndTests) {
   LL_UInt128 base({0xffff00000000ffff, 0xffffffff00000000});
   LL_UInt128 val128({0xf0f0f0f00f0f0f0f, 0xff00ff0000ff00ff});
   uint64_t val64 = 0xf0f0f0f00f0f0f0f;
   int val32 = 0x0f0f0f0f;
   LL_UInt128 result128({0xf0f0000000000f0f, 0xff00ff0000000000});
   LL_UInt128 result64(0xf0f0000000000f0f);
   LL_UInt128 result32(0x00000f0f);
   EXPECT_EQ((base & val128), result128);
   EXPECT_EQ((base & val64), result64);
   EXPECT_EQ((base & val32), result32);
 }
 
-TEST(LlvmLibcUInt128ClassTest, OrTests) {
+TEST(LlvmLibcUIntClassTest, OrTests) {
   LL_UInt128 base({0xffff00000000ffff, 0xffffffff00000000});
   LL_UInt128 val128({0xf0f0f0f00f0f0f0f, 0xff00ff0000ff00ff});
   uint64_t val64 = 0xf0f0f0f00f0f0f0f;
   int val32 = 0x0f0f0f0f;
   LL_UInt128 result128({0xfffff0f00f0fffff, 0xffffffff00ff00ff});
   LL_UInt128 result64({0xfffff0f00f0fffff, 0xffffffff00000000});
   LL_UInt128 result32({0xffff00000f0fffff, 0xffffffff00000000});
   EXPECT_EQ((base | val128), result128);
   EXPECT_EQ((base | val64), result64);
   EXPECT_EQ((base | val32), result32);
 }
 
-TEST(LlvmLibcUInt128ClassTest, CompoundAssignments) {
+TEST(LlvmLibcUIntClassTest, CompoundAssignments) {
   LL_UInt128 x({0xffff00000000ffff, 0xffffffff00000000});
   LL_UInt128 b({0xf0f0f0f00f0f0f0f, 0xff00ff0000ff00ff});
 
   LL_UInt128 a = x;
   a |= b;
   LL_UInt128 or_result({0xfffff0f00f0fffff, 0xffffffff00ff00ff});
   EXPECT_EQ(a, or_result);
 
TEST(LlvmLibcUIntClassTest, CompoundAssignments) {
 
   a = LL_UInt128({0xf, 0});
   b = LL_UInt128({0x1111111111111111, 0x1111111111111111});
   LL_UInt128 mul_result({0xffffffffffffffff, 0xffffffffffffffff});
   a *= b;
   EXPECT_EQ(a, mul_result);
 }
 
-TEST(LlvmLibcUInt128ClassTest, UnaryPredecrement) {
+TEST(LlvmLibcUIntClassTest, UnaryPredecrement) {
   LL_UInt128 a = LL_UInt128({0x1111111111111111, 0x1111111111111111});
   ++a;
   EXPECT_EQ(a, LL_UInt128({0x1111111111111112, 0x1111111111111111}));
 
   a = LL_UInt128({0xffffffffffffffff, 0x0});
   ++a;
   EXPECT_EQ(a, LL_UInt128({0x0, 0x1}));
 
   a = LL_UInt128({0xffffffffffffffff, 0xffffffffffffffff});
   ++a;
   EXPECT_EQ(a, LL_UInt128({0x0, 0x0}));
 }
 
-TEST(LlvmLibcUInt128ClassTest, EqualsTests) {
+TEST(LlvmLibcUIntClassTest, EqualsTests) {
   LL_UInt128 a1({0xffffffff00000000, 0xffff00000000ffff});
   LL_UInt128 a2({0xffffffff00000000, 0xffff00000000ffff});
   LL_UInt128 b({0xff00ff0000ff00ff, 0xf0f0f0f00f0f0f0f});
   LL_UInt128 a_reversed({0xffff00000000ffff, 0xffffffff00000000});
   LL_UInt128 a_upper(0xffff00000000ffff);
   LL_UInt128 a_lower(0xffffffff00000000);
   ASSERT_TRUE(a1 == a1);
   ASSERT_TRUE(a1 == a2);
   ASSERT_FALSE(a1 == b);
   ASSERT_FALSE(a1 == a_reversed);
   ASSERT_FALSE(a1 == a_lower);
   ASSERT_FALSE(a1 == a_upper);
   ASSERT_TRUE(a_lower != a_upper);
 }
 
-TEST(LlvmLibcUInt128ClassTest, ComparisonTests) {
+TEST(LlvmLibcUIntClassTest, ComparisonTests) {
   LL_UInt128 a({0xffffffff00000000, 0xffff00000000ffff});
   LL_UInt128 b({0xff00ff0000ff00ff, 0xf0f0f0f00f0f0f0f});
   EXPECT_GT(a, b);
   EXPECT_GE(a, b);
   EXPECT_LT(b, a);
   EXPECT_LE(b, a);
 
   LL_UInt128 x(0xffffffff00000000);
   LL_UInt128 y(0x00000000ffffffff);
   EXPECT_GT(x, y);
   EXPECT_GE(x, y);
   EXPECT_LT(y, x);
   EXPECT_LE(y, x);
 
   EXPECT_LE(a, a);
   EXPECT_GE(a, a);
 }
+
+TEST(LlvmLibcUIntClassTest, FullMulTests) {
+  LL_UInt128 a({0xffffffffffffffffULL, 0xffffffffffffffffULL});
+  LL_UInt128 b({0xfedcba9876543210ULL, 0xfefdfcfbfaf9f8f7ULL});
+  LL_UInt256 r({0x0123456789abcdf0ULL, 0x0102030405060708ULL,
+                0xfedcba987654320fULL, 0xfefdfcfbfaf9f8f7ULL});
+  LL_UInt128 r_hi({0xfedcba987654320eULL, 0xfefdfcfbfaf9f8f7ULL});
+
+  EXPECT_EQ(a.ful_mul(b), r);
+  EXPECT_EQ(a.quick_mul_hi(b), r_hi);
+
+  LL_UInt192 c(
+      {0x7766554433221101ULL, 0xffeeddccbbaa9988ULL, 0x1f2f3f4f5f6f7f8fULL});
+  LL_UInt320 rr({0x8899aabbccddeeffULL, 0x0011223344556677ULL,
+                 0x583715f4d3b29171ULL, 0xffeeddccbbaa9988ULL,
+                 0x1f2f3f4f5f6f7f8fULL});
+
+  EXPECT_EQ(a.ful_mul(c), rr);
+  EXPECT_EQ(a.ful_mul(c), c.ful_mul(a));
+}
+
+#define TEST_QUICK_MUL_HI(Bits, Error)                                         \
+  do {                                                                         \
+    LL_UInt##Bits a = ~LL_UInt##Bits(0);                                       \
+    LL_UInt##Bits hi = a.quick_mul_hi(a);                                      \
+    LL_UInt##Bits trunc = static_cast<LL_UInt##Bits>(a.ful_mul(a) >> Bits);    \
+    uint64_t overflow = trunc.sub(hi);                                         \
+    EXPECT_EQ(overflow, uint64_t(0));                                          \
+    EXPECT_LE(uint64_t(trunc), uint64_t(Error));                               \
+  } while (0)
+
+TEST(LlvmLibcUIntClassTest, QuickMulHiTests) {
+  TEST_QUICK_MUL_HI(128, 1);
+  TEST_QUICK_MUL_HI(192, 2);
+  TEST_QUICK_MUL_HI(256, 3);
+  TEST_QUICK_MUL_HI(512, 7);
+}

libc/utils/UnitTest/LibcTest.cpp

template bool test<__llvm_libc::cpp::UInt<192>>(
     __llvm_libc::cpp::UInt<192> RHS, const char *LHSStr, const char *RHSStr,
     const char *File, unsigned long Line);
 
 template bool test<__llvm_libc::cpp::UInt<256>>(
     RunContext *Ctx, TestCondition Cond, __llvm_libc::cpp::UInt<256> LHS,
     __llvm_libc::cpp::UInt<256> RHS, const char *LHSStr, const char *RHSStr,
     const char *File, unsigned long Line);
 
+template bool test<__llvm_libc::cpp::UInt<320>>(
+    RunContext *Ctx, TestCondition Cond, __llvm_libc::cpp::UInt<320> LHS,
+    __llvm_libc::cpp::UInt<320> RHS, const char *LHSStr, const char *RHSStr,
+    const char *File, unsigned long Line);
+
 template bool test<__llvm_libc::cpp::string_view>(
     RunContext *Ctx, TestCondition Cond, __llvm_libc::cpp::string_view LHS,
     __llvm_libc::cpp::string_view RHS, const char *LHSStr, const char *RHSStr,
     const char *File, unsigned long Line);
 
 } // namespace internal
 
 bool Test::testStrEq(const char *LHS, const char *RHS, const char *LHSStr,