sys/frac: Integer fractions

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