sys/div: Add support for big 64 bit numbers

Makefile.dep

@@ -535,6 +535,7 @@ endif
 
 ifneq (,$(filter xtimer,$(USEMODULE)))
     FEATURES_REQUIRED += periph_timer
+    USEMODULE += div
 endif
 
 ifneq (,$(filter saul_reg,$(USEMODULE)))

sys/div/Makefile

@@ -0,0 +1 @@
+include $(RIOTBASE)/Makefile.base

sys/div/div.c

@@ -0,0 +1,46 @@
+/**
+ * Copyright (C) 2016 Eistec AB
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser
+ * General Public License v2.1. See the file LICENSE in the top level
+ * directory for more details.
+ *
+ * @ingroup sys_util
+ * @{
+ * @file
+ * @brief    Integer division function implementations
+ *
+ * @author   Joakim Nohlgård <joakim.nohlgard@eistec.se>
+ *
+ * @}
+ */
+
+#include <stdint.h>
+#include <stdio.h>
+#include "div.h"
+
+uint64_t _div_mulhi64(const uint64_t a, const uint64_t b)
+{
+    /* Handle overflow explicit because we don't have 128 bit integers on
+     * our platforms. */
+    const uint32_t a_lo = (const uint32_t)a;
+    const uint32_t a_hi = (const uint32_t)(a >> 32);
+    const uint32_t b_lo = (const uint32_t)b;
+    const uint32_t b_hi = (const uint32_t)(b >> 32);
+
+    const uint64_t a_x_b_mid = (const uint64_t)a_hi * b_lo;
+    const uint64_t b_x_a_mid = (const uint64_t)b_hi * a_lo;
+    const uint64_t a_x_b_lo =  (const uint64_t)a_lo * b_lo;
+    const uint64_t a_x_b_hi =  (const uint64_t)a_hi * b_hi;
+
+    /* We may get up to 2 carry bits from the lower part of the multiplication */
+    const uint32_t carry_bits = ((uint64_t)(uint32_t)a_x_b_mid +
+                          (uint64_t)(uint32_t)b_x_a_mid +
+                          (a_x_b_lo >> 32) ) >> 32;
+
+    const uint64_t multhi = a_x_b_hi +
+                      (a_x_b_mid >> 32) + (b_x_a_mid >> 32) +
+                      carry_bits;
+
+    return multhi;
+}

sys/include/div.h

@@ -8,7 +8,7 @@
  */
 
 /**
- * @brief     Provides integer division functions
+ * @brief     Integer division functions
  *
  * This header provides some integer division functions that can be used
  * to prevent linking in compiler-generated ones, which are often larger.
@@ -31,35 +31,58 @@ extern "C" {
 #endif
 
 /**
- * @brief Integer divide val by 15625
+ * @brief Approximation of (2**l)/d for d=15625, l=12, 32 bits
+ */
+#define DIV_H_INV_15625_32    0x431bde83ul
+
+/**
+ * @brief Approximation of (2**l)/d for d=15625, l=12, 64 bits
+ */
+#define DIV_H_INV_15625_64    0x431bde82d7b634dbull
+
+/**
+ * @brief Required shifts for division by 15625, l above
+ */
+#define DIV_H_INV_15625_SHIFT 12
+
+/**
+ * @internal
+ * @brief Multiply two 64 bit integers into a 128 bit integer and return the upper half.
  *
- * @pre val <= 16383999997
+ * The implementation only uses 64 bit integers internally, no __int128 support
+ * is necessary.
+ *
+ * @see http://stackoverflow.com/questions/28868367/getting-the-high-part-of-64-bit-integer-multiplication
+
+ * @param[in]   a     operand a
+ * @param[in]   b     operand b
+ * @return      (((uint128_t)a * b) >> 64)
+ */
+uint64_t _div_mulhi64(const uint64_t a, const uint64_t b);
+
+/**
+ * @brief Integer divide val by 15625, 64 bit version
  *
  * @param[in]   val     dividend
  * @return      (val / 15625)
  */
 static inline uint64_t div_u64_by_15625(uint64_t val)
 {
-    /* a higher value would overflow 2^64 in the multiplication that follows */
-    assert(val <= 16383999997ull);
-
-    return (val * 0x431bde83UL) >> (12 + 32);
+    if (val > 16383999997ull) {
+        return (_div_mulhi64(DIV_H_INV_15625_64, val) >> DIV_H_INV_15625_SHIFT);
+    }
+    return (val * DIV_H_INV_15625_32) >> (DIV_H_INV_15625_SHIFT + 32);
 }
 
 /**
  * @brief Integer divide val by 1000000
  *
- * @pre val <= 1048575999808
- *
  * @param[in]   val     dividend
  * @return      (val / 1000000)
  */
 static inline uint64_t div_u64_by_1000000(uint64_t val)
 {
-    /* a higher value would overflow 2^64 in the multiplication that follows */
-    assert(val <= 1048575999808ull);
-
-    return div_u64_by_15625(val>>6);
+    return div_u64_by_15625(val) >> 6;
 }
 
 /**
@@ -68,15 +91,17 @@ static inline uint64_t div_u64_by_1000000(uint64_t val)
  * This is used to quantize a 1MHz value to the closest 32768Hz value,
  * e.g., for timers.
  *
- * The algorithm actually multiplies by 512 first, then divides by 15625,
- * keeping the result closer to a floored floating point division.
+ * The algorithm uses the modular multiplicative inverse of 15625 to use only
+ * multiplication and bit shifts to perform the division.
+ *
+ * The result will be equal to the mathematical expression: floor((val * 512) / 15625)
  *
  * @param[in]   val     dividend
  * @return      (val / (15625/512))
  */
 static inline uint32_t div_u32_by_15625div512(uint32_t val)
 {
-    return ((uint64_t)(val) * 0x431bde83ul) >> (12 + 32 - 9);
+    return ((uint64_t)(val) * DIV_H_INV_15625_32) >> (DIV_H_INV_15625_SHIFT + 32 - 9);
 }
 
 /**
@@ -85,23 +110,21 @@ static inline uint32_t div_u32_by_15625div512(uint32_t val)
  * This is used to quantize a 1MHz value to the closest 32768Hz value,
  * e.g., for timers.
  *
- * The algorithm actually multiplies by 512 first, then divides by 15625,
- * keeping the result closer to a floored floating point division.
- *
- * @pre val <= 16383999997
- *
  * @param[in]   val     dividend
  * @return      (val / (15625/512))
  */
 static inline uint64_t div_u64_by_15625div512(uint64_t val)
 {
-    /* a higher value would overflow 2^64 in the multiplication that follows */
-    assert(val <= 16383999997ull);
     /*
      * This saves around 1400 bytes of ROM on Cortex-M platforms (both ARMv6 and
      * ARMv7) from avoiding linking against __aeabi_uldivmod and related helpers
      */
-    return ((uint64_t)(val) * 0x431bde83ul) >> (12 + 32 - 9);
+    if (val > 16383999997ull) {
+        /* this would overflow 2^64 in the multiplication that follows, need to
+         * use the long version */
+        return (_div_mulhi64(DIV_H_INV_15625_64, val) >> (DIV_H_INV_15625_SHIFT - 9));
+    }
+    return (val * DIV_H_INV_15625_32) >> (DIV_H_INV_15625_SHIFT + 32 - 9);
 }
 
 /**

tests/unittests/tests-div/Makefile.include

@@ -0,0 +1 @@
+USEMODULE += div

tests/unittests/tests-div/tests-div.c

@@ -15,7 +15,7 @@
 #define ENABLE_DEBUG (0)
 #include "debug.h"
 
-static uint32_t u32_test_values[] = {
+static const uint32_t u32_test_values[] = {
     0,
     1,
     10,
@@ -28,26 +28,57 @@ static uint32_t u32_test_values[] = {
     0xffffffff,
 };
 
-static uint64_t u64_test_values[] = {
-    16383999997ull,
+static const uint64_t u64_test_values[] = {
     11111111111ull,
     0xffffffffull+1,
+    16383999997ull,
+    16383999998ull,
+    16383999999ull,
+    16384000000ull,
+    1048575999807ull,
+    1048575999808ull,
+    1048575999809ull,
+    0xffffffffffffeeull,
+    0x777777777777777ull,
+    0x1111111111111111ull,
+    0xffffffffffffffffull,
+    0x8000000000000000ull,
 };
 
-#define N_U32_VALS (sizeof(u32_test_values)/sizeof(uint32_t))
-#define N_U64_VALS (sizeof(u64_test_values)/sizeof(uint64_t))
+/* These expected values for test_div_u64_by_15625div512
+ * were computed from the expression (u64_test_values * 512) / 15625 using 128
+ * bit integers. */
+static const uint64_t u64_15625_512_expected_values[] = {
+    364088888,
+    140737488,
+    536870911,
+    536870911,
+    536870911,
+    536870912,
+    34359738361,
+    34359738361,
+    34359738361,
+    2361183241434822,
+    17630168202713342,
+    40297527320487639,
+    604462909807314587,
+    302231454903657293,
+};
+
+#define N_U32_VALS (sizeof(u32_test_values)/sizeof(u32_test_values[0]))
+#define N_U64_VALS (sizeof(u64_test_values)/sizeof(u64_test_values[0]))
 
 static void test_div_u64_by_15625(void)
 {
     for (unsigned i = 0; i < N_U32_VALS; i++) {
-        DEBUG("Dividing %"PRIu32" by 15625...\n", u32_test_values[i]);
+        DEBUG("Dividing %12"PRIu32" by 15625...\n", u32_test_values[i]);
         TEST_ASSERT_EQUAL_INT(
-            u32_test_values[i] / 15625,
+            (uint64_t)u32_test_values[i] / 15625,
             div_u64_by_15625(u32_test_values[i]));
     }
 
     for (unsigned i = 0; i < N_U64_VALS; i++) {
-        DEBUG("Dividing %"PRIu64" by 15625...\n", u64_test_values[i]);
+        DEBUG("Dividing %12"PRIu64" by 15625...\n", u64_test_values[i]);
         TEST_ASSERT_EQUAL_INT(
             (uint64_t)u64_test_values[i] / 15625,
             div_u64_by_15625(u64_test_values[i]));
@@ -69,7 +100,7 @@ static void test_div_u64_by_1000000(void)
     for (unsigned i = 0; i < N_U32_VALS; i++) {
         DEBUG("Dividing %"PRIu32" by 1000000...\n", u32_test_values[i]);
         TEST_ASSERT_EQUAL_INT(
-            u32_test_values[i] / 1000000lu,
+            (uint64_t)u32_test_values[i] / 1000000lu,
             div_u64_by_1000000(u32_test_values[i]));
     }
 
@@ -93,7 +124,7 @@ static void test_div_u64_by_15625div512(void)
     for (unsigned i = 0; i < N_U64_VALS; i++) {
         DEBUG("Dividing %"PRIu64" by (15625/512)...\n", u64_test_values[i]);
         TEST_ASSERT_EQUAL_INT(
-            (uint64_t)u64_test_values[i] * 512lu / 15625,
+            u64_15625_512_expected_values[i],
             div_u64_by_15625div512(u64_test_values[i]));
     }
 }