Improved fix16_exp() (half runtime, more accuracy), added fix16_log(), made unittests for both.

4 files changed, 217 insertions, 24 deletions

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 <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);

 }

+

+

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 <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;
+}
\ No newline at end of file