Improved fix16_exp() (half runtime, more accuracy), added fix16_log(), made unittests for both.
This commit is contained in:
Petteri.Aimonen
parent
eb1cb67561
commit
e929442f71
|
|
@ -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"
|
||||
|
|
|
|||
|
|
@ -1,23 +1,19 @@
|
|||
#include "fix16.h"
|
||||
#include "int64.h"
|
||||
|
||||
|
||||
#include <stdbool.h>
|
||||
|
||||
#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);
|
||||
}
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -0,0 +1,124 @@
|
|||
#include <fix16.h>
|
||||
#include <stdio.h>
|
||||
#include <math.h>
|
||||
#include <stdbool.h>
|
||||
#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;
|
||||
}
|
||||
Loading…
Reference in New Issue