diff --git a/libfixmath/fix16.h b/libfixmath/fix16.h index a2a49b1..31d7e00 100644 --- a/libfixmath/fix16.h +++ b/libfixmath/fix16.h @@ -159,6 +159,10 @@ extern fix16_t fix16_sqrt(fix16_t inValue) FIXMATH_FUNC_ATTRS; */ extern fix16_t fix16_exp(fix16_t inValue) FIXMATH_FUNC_ATTRS; +/*! Returns the natural logarithm of the given fix16_t. + */ +extern fix16_t fix16_log(fix16_t inValue) FIXMATH_FUNC_ATTRS; + #ifdef __cplusplus } #include "fix16.hpp" diff --git a/libfixmath/fix16_exp.c b/libfixmath/fix16_exp.c index e7e92e4..739ce04 100644 --- a/libfixmath/fix16_exp.c +++ b/libfixmath/fix16_exp.c @@ -1,23 +1,19 @@ #include "fix16.h" -#include "int64.h" - - +#include #ifndef FIXMATH_NO_CACHE static fix16_t _fix16_exp_cache_index[4096] = { 0 }; static fix16_t _fix16_exp_cache_value[4096] = { 0 }; #endif - - fix16_t fix16_exp(fix16_t inValue) { if(inValue == 0) return fix16_one; if(inValue == fix16_one) return fix16_e; - if(inValue > 681391) + if(inValue >= 681391) return fix16_max; - if(inValue < -726817) + if(inValue <= -772243) return 0; #ifndef FIXMATH_NO_CACHE @@ -26,23 +22,84 @@ fix16_t fix16_exp(fix16_t inValue) { if(_fix16_exp_cache_index[tempIndex] == inValue) return _fix16_exp_cache_value[tempIndex]; #endif - - int64_t tempOut = int64_add(int64_from_int32(fix16_one), int64_from_int32(inValue)); - int64_t tempValue = int64_from_int32(inValue); - uint32_t i, n; - for(i = 3, n = 2; i < 13; n *= i, i++) { - tempValue = int64_mul_i64_i32(tempValue, inValue); - #ifndef FIXMATH_NO_ROUNDING - tempValue = int64_add(tempValue, int64_from_int32(fix16_one >> 1)); - #endif - tempValue = int64_shift(tempValue, -16); - tempOut = int64_add(tempOut, int64_div_i64_i32(tempValue, n)); - } - + + /* The algorithm is based on the power series for exp(x): + * http://en.wikipedia.org/wiki/Exponential_function#Formal_definition + * + * From term n, we get term n+1 by multiplying with x/n. + * When the sum term drops to zero, we can stop summing. + */ + + // The power-series converges much faster on positive values + // and exp(-x) = 1/exp(x). + bool neg = (inValue < 0); + if (neg) + inValue = -inValue; + + fix16_t result = inValue + fix16_one; + fix16_t term = inValue; + int i; + + for (i = 2; i < 30; i++) + { + term = fix16_mul(term, fix16_div(inValue, fix16_from_int(i))); + result += term; + + if (term < 500 && (i > 15 || term < 20)) + break; + } + + if (neg) + result = fix16_div(fix16_one, result); + #ifndef FIXMATH_NO_CACHE _fix16_exp_cache_index[tempIndex] = inValue; - _fix16_exp_cache_value[tempIndex] = int64_lo(tempOut); + _fix16_exp_cache_value[tempIndex] = result; #endif - return int64_lo(tempOut); + return result; } + +fix16_t fix16_log(fix16_t inValue) +{ + fix16_t guess = fix16_from_int(2); + fix16_t delta; + int scaling = 0; + int count = 0; + + if (inValue <= 0) + return fix16_min; + + // Bring the value to the most accurate range (1 < x < 100) + const fix16_t e_to_fourth = 3578144; + while (inValue > fix16_from_int(100)) + { + inValue = fix16_div(inValue, e_to_fourth); + scaling += 4; + } + + while (inValue < fix16_from_int(1)) + { + inValue = fix16_mul(inValue, e_to_fourth); + scaling -= 4; + } + + do + { + // Solving e(x) = y using Newton's method + // f(x) = e(x) - y + // f'(x) = e(x) + fix16_t e = fix16_exp(guess); + delta = fix16_div(inValue - e, e); + + // It's unlikely that logarithm is very large, so avoid overshooting. + if (delta > fix16_from_int(3)) + delta = fix16_from_int(3); + + guess += delta; + } while (count++ < 10 && (delta > 1 || delta < -1)); + + return guess + fix16_from_int(scaling); +} + + diff --git a/unittests/Makefile b/unittests/Makefile index 80df66b..cb5dbb4 100644 --- a/unittests/Makefile +++ b/unittests/Makefile @@ -6,9 +6,9 @@ CFLAGS = -g -O0 -I../libfixmath -Wall -Wextra -Werror # The files required for tests FIX16_SRC = ../libfixmath/fix16.c ../libfixmath/fix16_sqrt.c \ - ../libfixmath/fix16.h + ../libfixmath/fix16_exp.c ../libfixmath/fix16.h -all: run_fix16_unittests +all: run_fix16_unittests run_fix16_exp_unittests clean: rm -f fix16_unittests_???? @@ -47,3 +47,11 @@ fix16_unittests_nn08: DEFINES=-DFIXMATH_NO_OVERFLOW -DFIXMATH_NO_ROUNDING -DFIXM fix16_unittests_% : fix16_unittests.c $(FIX16_SRC) $(CC) $(CFLAGS) $(DEFINES) -o $@ $^ -lm + + +# Tests for the exponential function, run only in default config +run_fix16_exp_unittests: fix16_exp_unittests + ./fix16_exp_unittests > /dev/null + +fix16_exp_unittests: fix16_exp_unittests.c $(FIX16_SRC) + $(CC) $(CFLAGS) $(DEFINES) -o $@ $^ -lm diff --git a/unittests/fix16_exp_unittests.c b/unittests/fix16_exp_unittests.c new file mode 100644 index 0000000..a033a05 --- /dev/null +++ b/unittests/fix16_exp_unittests.c @@ -0,0 +1,124 @@ +#include +#include +#include +#include +#include "unittests.h" + +#define delta(a,b) (((a)>=(b)) ? (a)-(b) : (b)-(a)) + +int main() +{ + int status = 0; + { + COMMENT("Testing fix16_exp() corner cases"); + TEST(fix16_exp(0) == fix16_one); + TEST(fix16_exp(fix16_min) == 0); + TEST(fix16_exp(fix16_max) == fix16_max); + } + + { + COMMENT("Testing fix16_exp() accuracy over -11..4"); + + fix16_t max_delta = -1; + fix16_t worst = 0; + fix16_t sum = 0; + int count = 0; + fix16_t a; + + for (a = fix16_from_dbl(-11.0); a < fix16_from_dbl(4.0); a += 31) + { + fix16_t result = fix16_exp(a); + fix16_t resultf = fix16_from_dbl(exp(fix16_to_dbl(a))); + + fix16_t d = delta(result, resultf); + if (d > max_delta) + { + max_delta = d; + worst = a; + } + + sum += d; + count++; + } + + printf("Worst delta %d with input %d\n", max_delta, worst); + printf("Average delta %0.2f\n", (float)sum / count); + + TEST(max_delta < 200); + } + + { + COMMENT("Testing fix16_exp() accuracy over full range"); + + float max_delta = -1; + fix16_t worst = 0; + float sum = 0; + int count = 0; + fix16_t a; + + // Test the whole range of results 0..32768 with a bit less samples + for (a = -772243; a < 681391; a += 113) + { + fix16_t result = fix16_exp(a); + fix16_t resultf = fix16_from_dbl(exp(fix16_to_dbl(a))); + + fix16_t d1 = delta(result, resultf); + + if (d1 > 0) d1--; // Forgive +-1 for the fix16_t inaccuracy + + float d = (float)d1 / resultf * 100; + + if (resultf < 1000) continue; // Percentages can explode when result is almost 0. + + if (d > max_delta) + { + max_delta = d; + worst = a; + } + + sum += d; + count++; + } + + printf("Worst delta %0.4f%% with input %d\n", max_delta, worst); + printf("Average delta %0.4f%%\n", sum / count); + + TEST(max_delta < 1); + } + + { + COMMENT("Testing fix16_log() accuracy over full range"); + + fix16_t max_delta = -1; + fix16_t worst = 0; + fix16_t sum = 0; + int count = 0; + fix16_t a; + + for (a = 100; a > 0 && a < fix16_max - 7561; a += 7561) + { + fix16_t result = fix16_log(a); + fix16_t resultf = fix16_from_dbl(log(fix16_to_dbl(a))); + + fix16_t d = delta(result, resultf); + if (d > max_delta) + { + max_delta = d; + worst = a; + } + + sum += d; + count++; + } + + printf("Worst delta %d with input %d\n", max_delta, worst); + printf("Average delta %0.2f\n", (float)sum / count); + + TEST(max_delta < 20); + } + + if (status != 0) + fprintf(stdout, "\n\nSome tests FAILED!\n"); + + return status; +} \ No newline at end of file