Add 128-bit multiplication and division operators

We already had addition and subtraction; this adds multiplication,
division, and modulus.

Makefile.am

@@ -263,7 +263,10 @@ u128_tests = \
     tests/u128-tests-shl.c.in \
     tests/u128-tests-shr.c.in \
     tests/u128-tests-add.c.in \
-    tests/u128-tests-sub.c.in
+    tests/u128-tests-sub.c.in \
+    tests/u128-tests-mul.c.in \
+    tests/u128-tests-div.c.in \
+    tests/u128-tests-mod.c.in
 
 CLEANFILES += $(u128_tests)
 
@@ -287,6 +290,12 @@ tests/u128-tests-add.c.in: tests/create-u128-test-cases.py
 	$(AM_V_GEN) $(PYTHON) $< add $@
 tests/u128-tests-sub.c.in: tests/create-u128-test-cases.py
 	$(AM_V_GEN) $(PYTHON) $< sub $@
+tests/u128-tests-mul.c.in: tests/create-u128-test-cases.py
+	$(AM_V_GEN) $(PYTHON) $< mul $@
+tests/u128-tests-div.c.in: tests/create-u128-test-cases.py
+	$(AM_V_GEN) $(PYTHON) $< div $@
+tests/u128-tests-mod.c.in: tests/create-u128-test-cases.py
+	$(AM_V_GEN) $(PYTHON) $< mod $@
 
 tests/test-u128.c: $(u128_tests)
 

include/libcork/core/u128.h

@@ -10,6 +10,7 @@
 #ifndef LIBCORK_CORE_U128_H
 #define LIBCORK_CORE_U128_H
 
+#include <stdlib.h>
 
 #include <libcork/config.h>
 #include <libcork/core/api.h>
@@ -268,6 +269,138 @@ cork_u128_sub(cork_u128 a, cork_u128 b)
     return result;
 }
 
+CORK_INLINE
+cork_u128
+cork_u128_mul(cork_u128 lhs, cork_u128 rhs)
+{
+    cork_u128  result;
+#if CORK_U128_HAVE_U128
+    result._.u128 = lhs._.u128 * rhs._.u128;
+#else
+    /* Multiplication table:  (each letter is 32 bits)
+     *
+     *         AaBb
+     *       * CcDd
+     *       ------
+     *           Ee = b * d
+     *          Ff  = B * d
+     *         Gg   = a * d
+     *        Hh    = A * d
+     *          Ii  = b * D
+     *         Jj   = B * D
+     *        Kk    = a * D
+     *       Ll     = A * D
+     *         Mm   = b * c
+     *        Nn    = B * c
+     *       Oo     = a * c         x = e
+     *      Pp      = A * c         X = E + f + i                 + [carry of x]
+     *        Qq    = b * C         w = F + g + I + j + m         + [carry of X]
+     *       Rr     = B * C         W = G + h + J + k + M + n + q + [carry of w]
+     *      Ss      = a * C         v = H + K + l + N + o + Q + r + [carry of W]
+     *     Tt       = A * C         V = L + O + p + R + s         + [carry of v]
+     *     --------                 u = P + S + t                 + [carry of V]
+     *     UuVvWwXx                 U = T                         + [carry of U]
+     */
+#endif
+    uint64_t A = lhs._.u32[3];
+    uint64_t a = lhs._.u32[2];
+    uint64_t B = lhs._.u32[1];
+    uint64_t b = lhs._.u32[0];
+    uint64_t C = rhs._.u32[3];
+    uint64_t c = rhs._.u32[2];
+    uint64_t D = rhs._.u32[1];
+    uint64_t d = rhs._.u32[0];
+    uint64_t Ee = b * d, E = Ee >> 32, e = Ee & 0xffffffff;
+    uint64_t Ff = B * d, F = Ff >> 32, f = Ff & 0xffffffff;
+    uint64_t Gg = a * d, G = Gg >> 32, g = Gg & 0xffffffff;
+    uint64_t Hh = A * d,               h = Hh & 0xffffffff;
+    uint64_t Ii = b * D, I = Ii >> 32, i = Ii & 0xffffffff;
+    uint64_t Jj = B * D, J = Jj >> 32, j = Jj & 0xffffffff;
+    uint64_t Kk = a * D,               k = Kk & 0xffffffff;
+    uint64_t Mm = b * c, M = Mm >> 32, m = Mm & 0xffffffff;
+    uint64_t Nn = B * c,               n = Nn & 0xffffffff;
+    uint64_t Qq = b * C,               q = Qq & 0xffffffff;
+    uint64_t x = e;
+    uint64_t X = E + f + i                 + (x >> 32);
+    uint64_t w = F + g + I + j + m         + (X >> 32);
+    uint64_t W = G + h + J + k + M + n + q + (w >> 32);
+    result._.u32[3] = W;
+    result._.u32[2] = w;
+    result._.u32[1] = X;
+    result._.u32[0] = x;
+    return result;
+}
+
+struct cork_u128_divmod {
+    cork_u128  div;
+    cork_u128  mod;
+};
+
+CORK_INLINE
+struct cork_u128_divmod
+cork_u128_divmod(cork_u128 a, cork_u128 b)
+{
+    struct cork_u128_divmod  result;
+#if CORK_U128_HAVE_U128
+    result.div._.u128 = a._.u128 / b._.u128;
+    result.mod._.u128 = a._.u128 % b._.u128;
+#else
+    cork_u128  divisor;
+    cork_u128  multiple;
+
+    if (cork_u128_eq(b, cork_u128_zero())) {
+        // Division by 0
+        abort();
+    }
+
+    divisor = b;
+    result.mod = a;
+    result.div = cork_u128_zero();
+    multiple = cork_u128_from_64(0, 1);
+
+    while (cork_u128_lt(divisor, a)) {
+        divisor = cork_u128_shl(divisor, 1);
+        multiple = cork_u128_shr(multiple, 1);
+    }
+
+    do {
+        if (cork_u128_ge(result.mod, divisor)) {
+            result.mod = cork_u128_sub(result.mod, divisor);
+            result.div = cork_u128_add(result.div, multiple);
+        }
+        divisor = cork_u128_shr(divisor, 1);
+        multiple = cork_u128_shr(multiple, 1);
+    } while (!cork_u128_eq(multiple, cork_u128_zero()));
+#endif
+    return result;
+}
+
+CORK_INLINE
+cork_u128
+cork_u128_div(cork_u128 a, cork_u128 b)
+{
+#if CORK_U128_HAVE_U128
+    cork_u128  result;
+    result._.u128 = a._.u128 / b._.u128;
+    return result;
+#else
+    return cork_u128_divmod(a, b).div;
+#endif
+}
+
+CORK_INLINE
+cork_u128
+cork_u128_mod(cork_u128 a, cork_u128 b)
+{
+#if CORK_U128_HAVE_U128
+    cork_u128  result;
+    result._.u128 = a._.u128 % b._.u128;
+    return result;
+#else
+    return cork_u128_divmod(a, b).div;
+#endif
+}
+
 
 /* log10(x) = log2(x) / log2(10) ~= log2(x) / 3.322 */
 #define CORK_U128_DECIMAL_LENGTH  44  /* ~= 128 / 3 + 1 + 1 */

src/libcork/core/u128.c

@@ -126,3 +126,15 @@ cork_u128_add(cork_u128 a, cork_u128 b);
 
 cork_u128
 cork_u128_sub(cork_u128 a, cork_u128 b);
+
+cork_u128
+cork_u128_mul(cork_u128 lhs, cork_u128 rhs);
+
+struct cork_u128_divmod
+cork_u128_divmod(cork_u128 a, cork_u128 b);
+
+cork_u128
+cork_u128_div(cork_u128 a, cork_u128 b);
+
+cork_u128
+cork_u128_mod(cork_u128 a, cork_u128 b);

tests/CMakeLists.txt

@@ -77,6 +77,9 @@ make_u128_suite(shl)
 make_u128_suite(shr)
 make_u128_suite(add)
 make_u128_suite(sub)
+make_u128_suite(mul)
+make_u128_suite(div)
+make_u128_suite(mod)
 
 #-----------------------------------------------------------------------
 # Command-line tests

tests/create-u128-test-cases.py

@@ -124,6 +124,62 @@ def create_one_sub_test_case():
     print("},", sep="")
 
 
+def create_one_mul_test_case():
+    lhs = random_128()
+    rhs = random_128()
+    result = (lhs * rhs) % 2**128
+    print()
+    print("/*    ", dec_128(lhs), sep="")
+    print(" *  * ", dec_128(rhs), sep="")
+    print(" *  = ", dec_128(result), sep="")
+    print(" */")
+    print("{", sep="")
+    print("    ", hex_128(lhs), ",", sep="")
+    print("    ", hex_128(rhs), ",", sep="")
+    print("    ", hex_128(result), sep="")
+    print("},", sep="")
+
+
+def create_one_div_test_case():
+    lhs = random_128()
+    rhs = random_128()
+    if rhs == 0:
+        rhs = 1
+    if lhs < rhs:
+        lhs, rhs = rhs, lhs
+    result = (lhs // rhs) % 2**128
+    print()
+    print("/*    ", dec_128(lhs), sep="")
+    print(" *  / ", dec_128(rhs), sep="")
+    print(" *  = ", dec_128(result), sep="")
+    print(" */")
+    print("{", sep="")
+    print("    ", hex_128(lhs), ",", sep="")
+    print("    ", hex_128(rhs), ",", sep="")
+    print("    ", hex_128(result), sep="")
+    print("},", sep="")
+
+
+def create_one_mod_test_case():
+    lhs = random_128()
+    rhs = random_128()
+    if rhs == 0:
+        rhs = 1
+    if lhs < rhs:
+        lhs, rhs = rhs, lhs
+    result = (lhs % rhs) % 2**128
+    print()
+    print("/*    ", dec_128(lhs), sep="")
+    print(" *  % ", dec_128(rhs), sep="")
+    print(" *  = ", dec_128(result), sep="")
+    print(" */")
+    print("{", sep="")
+    print("    ", hex_128(lhs), ",", sep="")
+    print("    ", hex_128(rhs), ",", sep="")
+    print("    ", hex_128(result), sep="")
+    print("},", sep="")
+
+
 if len(sys.argv) == 1:
     print("Usage: create-u128-test-cases.py [operator]")
     sys.exit(1)
@@ -154,3 +210,9 @@ def create_one_sub_test_case():
         create_one_add_test_case()
     elif sys.argv[1] == "sub":
         create_one_sub_test_case()
+    elif sys.argv[1] == "mul":
+        create_one_mul_test_case()
+    elif sys.argv[1] == "div":
+        create_one_div_test_case()
+    elif sys.argv[1] == "mod":
+        create_one_mod_test_case()

tests/test-u128.c

@@ -304,6 +304,60 @@ START_TEST(test_u128_sub)
 END_TEST
 
 
+static const struct arithmetic_test MUL_TESTS[] = {
+    {0, 0, 0, 0, 0, 0},
+    {0, 1, 0, 1, 0, 1},
+    {0, 2, 0, 1, 0, 2},
+    {1, 0, 0, 1, 1, 0},
+    {1, 0, 0, 2, 2, 0},
+    {0, UINT64_C(0x8000000000000000), 0, 2, 1, 0},
+#include "u128-tests-mul.c.in"
+};
+
+START_TEST(test_u128_mul)
+{
+    DESCRIBE_TEST;
+    check_arithmetic_tests(cork_u128_mul, "*", MUL_TESTS);
+}
+END_TEST
+
+
+static const struct arithmetic_test DIV_TESTS[] = {
+    {0, 0, 0, 1, 0, 0},
+    {0, 1, 0, 1, 0, 1},
+    {0, 2, 0, 1, 0, 2},
+    {1, 0, 0, 1, 1, 0},
+    {2, 0, 0, 2, 1, 0},
+    {1, 0, 0, 2, 0, UINT64_C(0x8000000000000000)},
+#include "u128-tests-div.c.in"
+};
+
+START_TEST(test_u128_div)
+{
+    DESCRIBE_TEST;
+    check_arithmetic_tests(cork_u128_div, "/", DIV_TESTS);
+}
+END_TEST
+
+
+static const struct arithmetic_test MOD_TESTS[] = {
+    {0, 0, 0, 1, 0, 0},
+    {0, 1, 0, 1, 0, 0},
+    {0, 2, 0, 1, 0, 0},
+    {1, 0, 0, 1, 0, 0},
+    {0, 3, 0, 2, 0, 1},
+    {1, 3, 0, 2, 0, 1},
+#include "u128-tests-mod.c.in"
+};
+
+START_TEST(test_u128_mod)
+{
+    DESCRIBE_TEST;
+    check_arithmetic_tests(cork_u128_mod, "%", MOD_TESTS);
+}
+END_TEST
+
+
 struct comparison_test {
     uint64_t i0;
     uint64_t i1;
@@ -454,6 +508,9 @@ test_suite()
     tcase_add_test(tc_u128, test_u128_shr);
     tcase_add_test(tc_u128, test_u128_add);
     tcase_add_test(tc_u128, test_u128_sub);
+    tcase_add_test(tc_u128, test_u128_mul);
+    tcase_add_test(tc_u128, test_u128_div);
+    tcase_add_test(tc_u128, test_u128_mod);
     tcase_add_test(tc_u128, test_u128_eq);
     tcase_add_test(tc_u128, test_u128_ne);
     tcase_add_test(tc_u128, test_u128_lt);