From e4e860a20a6b7a881ec64bd384aed2a32c3024ff Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Joakim=20Nohlg=C3=A5rd?= Date: Tue, 4 Sep 2018 15:13:44 +0200 Subject: [PATCH] frac: Integer fraction scaling library --- sys/frac/Makefile | 1 + sys/frac/frac.c | 157 +++++++++++++++++++ sys/include/frac.h | 91 +++++++++++ tests/frac-config/Makefile | 9 ++ tests/frac-config/README.md | 17 +++ tests/frac-config/main.c | 54 +++++++ tests/unittests/tests-frac/Makefile | 1 + tests/unittests/tests-frac/Makefile.include | 1 + tests/unittests/tests-frac/tests-frac.c | 158 ++++++++++++++++++++ tests/unittests/tests-frac/tests-frac.h | 43 ++++++ 10 files changed, 532 insertions(+) create mode 100644 sys/frac/Makefile create mode 100644 sys/frac/frac.c create mode 100644 sys/include/frac.h create mode 100644 tests/frac-config/Makefile create mode 100644 tests/frac-config/README.md create mode 100644 tests/frac-config/main.c create mode 100644 tests/unittests/tests-frac/Makefile create mode 100644 tests/unittests/tests-frac/Makefile.include create mode 100644 tests/unittests/tests-frac/tests-frac.c create mode 100644 tests/unittests/tests-frac/tests-frac.h diff --git a/sys/frac/Makefile b/sys/frac/Makefile new file mode 100644 index 0000000000..48422e909a --- /dev/null +++ b/sys/frac/Makefile @@ -0,0 +1 @@ +include $(RIOTBASE)/Makefile.base diff --git a/sys/frac/frac.c b/sys/frac/frac.c new file mode 100644 index 0000000000..35c69855e6 --- /dev/null +++ b/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 + * + * @} + */ + +#include +#include +#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); +} diff --git a/sys/include/frac.h b/sys/include/frac.h new file mode 100644 index 0000000000..76738b32f4 --- /dev/null +++ b/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 + * @{ + */ + +#ifndef FRAC_H +#define FRAC_H + +#include + +#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 */ diff --git a/tests/frac-config/Makefile b/tests/frac-config/Makefile new file mode 100644 index 0000000000..1490983b30 --- /dev/null +++ b/tests/frac-config/Makefile @@ -0,0 +1,9 @@ +include ../Makefile.tests_common + +BOARD ?= native + +BOARD_WHITELIST += native + +USEMODULE += frac + +include $(RIOTBASE)/Makefile.include diff --git a/tests/frac-config/README.md b/tests/frac-config/README.md new file mode 100644 index 0000000000..ee41772e34 --- /dev/null +++ b/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); diff --git a/tests/frac-config/main.c b/tests/frac-config/main.c new file mode 100644 index 0000000000..9c60c95c15 --- /dev/null +++ b/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 + * + * @} + */ + +#include +#include +#include + +#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; +} diff --git a/tests/unittests/tests-frac/Makefile b/tests/unittests/tests-frac/Makefile new file mode 100644 index 0000000000..48422e909a --- /dev/null +++ b/tests/unittests/tests-frac/Makefile @@ -0,0 +1 @@ +include $(RIOTBASE)/Makefile.base diff --git a/tests/unittests/tests-frac/Makefile.include b/tests/unittests/tests-frac/Makefile.include new file mode 100644 index 0000000000..74f7ec33c5 --- /dev/null +++ b/tests/unittests/tests-frac/Makefile.include @@ -0,0 +1 @@ +USEMODULE += frac diff --git a/tests/unittests/tests-frac/tests-frac.c b/tests/unittests/tests-frac/tests-frac.c new file mode 100644 index 0000000000..525ffa10c7 --- /dev/null +++ b/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 +#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()); +} diff --git a/tests/unittests/tests-frac/tests-frac.h b/tests/unittests/tests-frac/tests-frac.h new file mode 100644 index 0000000000..ea4db48b5f --- /dev/null +++ b/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 + */ +#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 */ +/** @} */