frac: Integer fraction scaling library

sys/frac/Makefile

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

sys/frac/frac.c

@@ -0,0 +1,157 @@
+/**
+ * Copyright (C) 2018 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_frac
+ * @{
+ * @file
+ * @brief    Integer fraction function implementations
+ *
+ * @author   Joakim Nohlgård <joakim.nohlgard@eistec.se>
+ *
+ * @}
+ */
+
+#include <stdint.h>
+#include <stdio.h>
+#include "assert.h"
+#include "frac.h"
+
+#define ENABLE_DEBUG (0)
+#include "debug.h"
+
+/**
+ * @brief   compute greatest common divisor of @p u and @p v
+ *
+ * @param[in]   u    first operand
+ * @param[in]   v    second operand
+ *
+ * @return  Greatest common divisor of @p u and @p v
+ */
+static uint32_t gcd32(uint32_t u, uint32_t v)
+{
+    /* Source: https://en.wikipedia.org/wiki/Binary_GCD_algorithm#Iterative_version_in_C */
+    unsigned shift;
+
+    /* GCD(0,v) == v; GCD(u,0) == u, GCD(0,0) == 0 */
+    if (u == 0) {
+        return v;
+    }
+    if (v == 0) {
+        return u;
+    }
+
+    /* Let shift := log2 K, where K is the greatest power of 2
+     * dividing both u and v. */
+    for (shift = 0; ((u | v) & 1) == 0; ++shift) {
+        u >>= 1;
+        v >>= 1;
+    }
+
+    /* remove all factors of 2 in u */
+    while ((u & 1) == 0) {
+        u >>= 1;
+    }
+
+    /* From here on, u is always odd. */
+    do {
+        /* remove all factors of 2 in v -- they are not common */
+        /*   note: v is not zero, so while will terminate */
+        while ((v & 1) == 0) {
+            v >>= 1;
+        }
+
+        /* Now u and v are both odd. Swap if necessary so u <= v,
+         * then set v = v - u (which is even). */
+        if (u > v) {
+            /* Swap u and v */
+            uint32_t t = v;
+            v = u;
+            u = t;
+        }
+
+        v = v - u; /* Here v >= u */
+    } while (v != 0);
+
+    /* restore common factors of 2 */
+    return u << shift;
+}
+
+uint32_t frac_long_divide(uint32_t num, uint32_t den, int *prec, uint32_t *rem)
+{
+    /* Binary long division with adaptive number of fractional bits */
+    /* The result will be a Qx.y number where x is the number of bits in the
+     * integer part and y = 64 - x. Similar to floating point, except the result
+     * is unsigned, and we can only represent numbers in the range 2**-32..(2**32 - 1) */
+    assert(den); /* divide by zero */
+
+    uint32_t q = 0; /* Quotient */
+    uint64_t r = 0;  /* Remainder */
+    if (prec) {
+        *prec = 0;
+    }
+    if (num == 0) {
+        if (rem) {
+            *rem = 0;
+        }
+        return 0;
+    }
+    unsigned p = bitarithm_msb(num);
+    int i_bits = p + 1; /* Number of integer bits in the result */
+    uint32_t num_mask = (1ul << p);
+    for (unsigned k = 0; k < (64u + p); ++k) {
+        r <<= 1;
+        q <<= 1;
+        if (num & num_mask) {
+            r |= 1;
+        }
+        num_mask >>= 1;
+        if (r >= den) {
+            r -= den;
+            q |= 1;
+        }
+        if (q == 0) {
+            --i_bits;
+        }
+        if (q & (1ul << 31u)) {
+            /* result register is full */
+            break;
+        }
+        if ((r == 0) && (num == 0)) {
+            /* divides evenly */
+            break;
+        }
+    }
+    if (r > 0) {
+        ++q;
+    }
+    if (prec) {
+        *prec = i_bits;
+    }
+    if (rem) {
+        *rem = r;
+    }
+
+    return q;
+}
+
+void frac_init(frac_t *frac, uint32_t num, uint32_t den)
+{
+    DEBUG("frac_init32(%p, %" PRIu32 ", %" PRIu32 ")\n", (const void *)frac, num, den);
+    assert(den);
+    /* Reduce the fraction to shortest possible form by dividing by the greatest
+     * common divisor */
+    uint32_t gcd = gcd32(num, den);
+    /* Divide den and num by their greatest common divisor */
+    den /= gcd;
+    num /= gcd;
+    int prec = 0;
+    uint32_t rem = 0;
+    frac->frac = frac_long_divide(num, den, &prec, &rem);
+    frac->shift = (sizeof(frac->frac) * 8) - prec;
+    DEBUG("frac_init32: gcd = %" PRIu32 " num = %" PRIu32 " den = %" PRIu32 " frac = 0x%08" PRIx32 " shift = %02d, rem = 0x%08" PRIx32 "\n",
+          gcd, num, den, frac->frac, frac->shift, rem);
+}

sys/include/frac.h

@@ -0,0 +1,91 @@
+/*
+ * Copyright (C) 2018 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.
+ */
+
+/**
+ * @defgroup  sys_frac   Fractional integer operations
+ * @ingroup   sys
+ *
+ * This header provides some functions for scaling integers by fractions, while
+ * preserving as many bits as possible.
+ *
+ * The implementation requires that @ref frac_t is initialized properly, either
+ * by calling @ref frac_init, which will compute the algorithm parameters at
+ * runtime, or via a precomputed initializer.
+ *
+ * Precomputing the frac_t values can be done via the application found in
+ * `tests/frac-config` in the RIOT tree.
+ *
+ * ### Numeric precision
+ *
+ * The algorithm will under certain circumstances give an incorrectly rounded
+ * result, more precisely, the result may sometimes be rounded up instead of
+ * rounded down when the product in the numerator, @$p = x * num@$, would
+ * result in @$p >= 2^{31}@$. Fortunately, the relative error of this rounding
+ * mistake is small.
+ *
+ * This tradeoff is a design choice to make the algorithm faster.
+ *
+ * @see       Libdivide homepage: http://libdivide.com/
+ *
+ * @file
+ * @ingroup   sys
+ * @author    Joakim Nohlgård <joakim.nohlgard@eistec.se>
+ * @{
+ */
+
+#ifndef FRAC_H
+#define FRAC_H
+
+#include <stdint.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/**
+ * @brief   frac descriptor for fraction consisting of two 32 bit integers
+ */
+typedef struct {
+    uint32_t frac; /**< fraction */
+    uint8_t shift; /**< exponent */
+} frac_t;
+
+/**
+ * @brief   Initialize frac_t struct
+ *
+ * This function computes the mathematical parameters used by the frac algorithm.
+ *
+ * @note If @p num > @p den, the result from @ref frac_scale modulo 2**32.
+ *
+ * @pre @p den must not be 0
+ *
+ * @param[out]  frac    pointer to frac descriptor to initialize
+ * @param[in]   num     numerator
+ * @param[in]   den     denominator
+ */
+void frac_init(frac_t *frac, uint32_t num, uint32_t den);
+
+/**
+ * @brief   Scale a 32 bit integer by a 32/32 rational number
+ *
+ * @param[in]   frac  scaling fraction
+ * @param[in]   x     unscaled integer
+ *
+ * @return      (x * frac) % 2**32, avoiding truncation
+ */
+static inline uint32_t frac_scale(const frac_t *frac, uint32_t x)
+{
+    uint32_t scaled = ((uint64_t)frac->frac * x) >> frac->shift;
+    return scaled;
+}
+
+#ifdef __cplusplus
+}
+#endif
+/** @} */
+#endif /* FRAC_H */

tests/frac-config/Makefile

@@ -0,0 +1,9 @@
+include ../Makefile.tests_common
+
+BOARD ?= native
+
+BOARD_WHITELIST += native
+
+USEMODULE += frac
+
+include $(RIOTBASE)/Makefile.include

tests/frac-config/README.md

@@ -0,0 +1,17 @@
+This is a helper program to generate frac_t static initializers, to avoid the
+runtime overhead of calling frac_init during initialization.
+
+When using the precomputed values it is possible to declare const frac_t in ROM,
+which is not possible with runtime computation.
+
+Example static configuration:
+
+    static const frac_t myfrac = {
+        .num = 512,
+        .den = 15625,
+        .div = { .magic = 0x8637bd05af6c69b6ull, .more = 0x0d },
+    };
+
+you can use myfrac as usual:
+
+    scaled_number = frac_scale(&myfrac, unscaled_number);

tests/frac-config/main.c

@@ -0,0 +1,54 @@
+/*
+ * Copyright (C) 2018 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     tests
+ * @{
+ *
+ * @file
+ * @brief       Frac library static configuration helper
+ *
+ * @author      Joakim Nohlgård <joakim.nohlgard@eistec.se>
+ *
+ * @}
+ */
+
+#include <stdio.h>
+#include <stdint.h>
+#include <inttypes.h>
+
+#include "frac.h"
+
+int main(void)
+{
+    puts("frac library static configuration generator");
+    while (1) {
+        int res = 0;
+        uint32_t num = 0;
+        uint32_t den = 0;
+        do {
+            puts("Enter fraction numerator (multiplier):");
+            res = scanf("%" SCNu32, &num);
+        } while (res == 0);
+        do {
+            puts("Enter fraction denominator (divider):");
+            res = scanf("%" SCNu32, &den);
+        } while (res == 0);
+        if (den == 0) {
+            continue;
+        }
+        frac_t frac;
+        frac_init(&frac, num, den);
+
+        printf("Static initialization of frac_t for (%" PRIu32 " / %" PRIu32 "):\n",
+            num, den);
+        puts("(Copy and paste into your code)\n");
+        printf("{ .frac = 0x%" PRIx32 ", .shift = %u }\n\n", frac.frac, (unsigned)frac.shift);
+    }
+    return 0;
+}

tests/unittests/tests-frac/Makefile

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

tests/unittests/tests-frac/Makefile.include

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

tests/unittests/tests-frac/tests-frac.c

@@ -0,0 +1,158 @@
+/*
+ * Copyright (C) 2018 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.
+ */
+
+#include <string.h>
+#include "embUnit.h"
+#include "tests-frac.h"
+
+#include "kernel_defines.h"
+#include "frac.h"
+#include "div.h"
+
+#define ENABLE_DEBUG (0)
+#include "debug.h"
+
+static const uint32_t u32_fraction_operands[] = {
+    1ul,
+    2ul,
+    5ul,
+    10ul,
+    100ul,
+    1000ul,
+    1000000ul,
+    2000000ul,
+    4000000ul,
+    8000000ul,
+    16000000ul,
+    32000000ul,
+    641ul,
+    274177ul,
+    32768ul,
+    9600ul,
+    38400ul,
+    115200ul,
+    230400ul,
+    460800ul,
+    921600ul,
+    4096ul,
+    15625ul,
+    125ul,
+    1048576ul,
+    0x10000000ul,
+    0x1000000ul,
+    1000000000ul,
+    999999733ul,  /* <- prime */
+    512000000ul,
+    1024000000ul,
+    0x40000000ul,
+    0x80000000ul,
+    0xc0000000ul,
+    0xe0000000ul,
+    0xf0000000ul,
+    0xfffffff0ul,
+    0xfffffff8ul,
+    0xfffffffcul,
+    0xfffffffeul,
+    0xfffffffful,
+};
+
+static const uint32_t u32_test_values[] = {
+    0ul,
+    1ul,
+    10ul,
+    32ul,
+    15625ul,
+    15625ul*5,
+    (15625ul*5)+1,
+    0xfffful,
+    0xfffful<<10,
+    1234567890ul,
+    99999999ul,
+    1000000ul,
+    115200ul,
+    38400ul,
+    57600ul,
+    921600ul,
+    32768ul,
+    16000000ul,
+    15999999ul,
+    32767ul,
+    327679999ul,
+    100000000ul,
+    2100012683ul,            /* <- prime */
+    0x7ffffffful,
+    0x80000000ul,
+    0xc0000000ul,
+    0xe0000000ul,
+    0xf0000000ul,
+    0xfffffff0ul,
+    0xfffffff8ul,
+    0xfffffffcul,
+    0xfffffffeul,
+    0xfffffffful,
+};
+
+#define N_U32_OPERANDS ARRAY_SIZE(u32_fraction_operands)
+#define N_U32_VALS ARRAY_SIZE(u32_test_values)
+
+static void test_frac_scale32(void)
+{
+    for (unsigned k = 0; k < N_U32_OPERANDS; ++k) {
+        for (unsigned j = 0; j < N_U32_OPERANDS; ++j) {
+            uint32_t num = u32_fraction_operands[j];
+            uint32_t den = u32_fraction_operands[k];
+            frac_t frac;
+            frac_init(&frac, num, den);
+            for (unsigned i = 0; i < N_U32_VALS; i++) {
+                DEBUG("Scaling %" PRIu32 " by (%" PRIu32 " / %" PRIu32 "), ",
+                    u32_test_values[i], num, den);
+                /* intermediate result */
+                volatile uint64_t tmp = (uint64_t)u32_test_values[i] * num;
+                volatile uint64_t expected = tmp / (uint64_t)den;
+                if (expected > 0xfffffffful) {
+                    expected &= 0xfffffffful;
+                }
+                uint32_t actual = frac_scale(&frac, u32_test_values[i]);
+                if ((uint32_t)expected != actual) {
+                    int32_t diff = actual - expected;
+                    DEBUG("%" PRIu32 " * (%" PRIu32 " / %" PRIu32 ")"
+                        " tmp %" PRIu64 " expect %" PRIu32 ", actual %" PRIu32 ", diff = %" PRId32 " shift=%u\n",
+                        u32_test_values[i], num, den, tmp, (uint32_t)expected,
+                        actual, diff, frac.shift);
+
+                    /* The frac algorithm sacrifices accuracy for speed,
+                     * some large numbers will be incorrectly rounded,
+                     * resulting in small differences here.. */
+                    uint32_t max_error = frac_scale(&frac, 2);
+                    max_error = max_error ? max_error : 1;
+                    TEST_ASSERT_EQUAL_INT(1, diff >= 0);
+                    if ((uint32_t)diff <= max_error) {
+                        continue;
+                    }
+                }
+                TEST_ASSERT_EQUAL_INT((uint32_t)expected, (uint32_t)actual);
+            }
+        }
+    }
+}
+
+Test *tests_frac_tests(void)
+{
+    EMB_UNIT_TESTFIXTURES(fixtures) {
+        new_TestFixture(test_frac_scale32),
+    };
+
+    EMB_UNIT_TESTCALLER(frac_tests, NULL, NULL, fixtures);
+
+    return (Test *)&frac_tests;
+}
+
+void tests_frac(void)
+{
+    TESTS_RUN(tests_frac_tests());
+}

tests/unittests/tests-frac/tests-frac.h

@@ -0,0 +1,43 @@
+/*
+ * Copyright (C) 2018 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.
+ */
+
+/**
+ * @addtogroup  unittests
+ * @{
+ *
+ * @file
+ * @brief       Unittests for the ``frac`` header
+ *
+ * @author      Joakim Nohlgård <joakim.nohlgard@eistec.se>
+ */
+#ifndef TESTS_FRAC_H
+#define TESTS_FRAC_H
+#include "embUnit/embUnit.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/**
+*  @brief   The entry point of this test suite.
+*/
+void tests_frac(void);
+
+/**
+ * @brief   Generates tests for frac
+ *
+ * @return  embUnit tests if successful, NULL if not.
+ */
+Test *tests_frac_tests(void);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* TESTS_FRAC_H */
+/** @} */