Added a test to assess the quality of the approximate math functions.

.gitignore

@@ -53,6 +53,7 @@ tests/testSingle
 tests/testTimeIntegration
 tests/testSPHStep
 tests/testKernel
+tests/testMaths
 tests/testParser
 tests/parser_output.yml
 tests/test27cells.sh

src/approx_math.h

@@ -24,12 +24,16 @@
 /**
  * @brief Approximate version of expf(x) using a 4th order Taylor expansion
  *
- * The absolute error is of order 10^-6 for -0.2 < x < 0.2.
+ * The absolute error is smaller than 3 * 10^-6 for -0.2 < x < 0.2.
+ * The absolute error is smaller than 2 * 10^-7 for -0.1 < x < 0.1.
+
+ * The relative error is smaller than 1 * 10^-6 for -0.2 < x < 0.2.
+ * The relative error is smaller than 4 * 10^-8 for -0.1 < x < 0.1.
  *
  * @param x The number to take the exponential of.
  */
 __attribute__((always_inline)) INLINE static float approx_expf(float x) {
-  return 1.f + x * (1.f + x * (0.5f + x * (1.f / 6.0f + 1.f / 24.0f * x)));
+  return 1.f + x * (1.f + x * (0.5f + x * (((float)(1.0 / 6.0)) + ((float)(1.0 / 24.0)) * x)));
 }
 
 #endif /* SWIFT_APPROX_MATH_H */

tests/Makefile.am

@@ -21,16 +21,18 @@ AM_CFLAGS = -I$(top_srcdir)/src $(HDF5_CPPFLAGS) -DTIMER
 AM_LDFLAGS = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS)
 
 # List of programs and scripts to run in the test suite
-TESTS = testGreetings testReading.sh testSingle testPair.sh testPairPerturbed.sh \
+TESTS = testGreetings testMaths testReading.sh testSingle testPair.sh testPairPerturbed.sh \
 	test27cells.sh test27cellsPerturbed.sh testParser.sh testKernel testSPHStep
 
 # List of test programs to compile
-check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
+check_PROGRAMS = testGreetings testMaths testReading testSingle testTimeIntegration \
 		 testSPHStep testPair test27cells testParser testKernel testInteractions
 
 # Sources for the individual programs
 testGreetings_SOURCES = testGreetings.c
 
+testMaths_SOURCES = testMaths.c
+
 testReading_SOURCES = testReading.c
 
 testTimeIntegration_SOURCES = testTimeIntegration.c

tests/testMaths.c

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#include "../config.h"
+
+#include "approx_math.h"
+#include "vector.h"
+
+#include <stdio.h>
+#include <math.h>
+
+int main() {
+
+  const int numPoints = 60000;
+  
+  for (int i = 0; i < numPoints; ++i) {
+    
+    const float x =  0.6f * (i / (float) numPoints) - 0.3f;
+    const float exp_correct = expf(x);
+    const float exp_approx = approx_expf(x);
+
+    const float abs = fabs(exp_correct - exp_approx);
+    const float rel = 0.5f * fabs(exp_correct - exp_approx) / fabs(exp_correct + exp_approx);
+    
+    printf("%2d: x= %f exp(x)= %e approx_exp(x)=%e abs=%e rel=%e\n", i, x, exp_correct,
+	   exp_approx, abs, rel);
+
+    if (abs > 3e-6 && fabsf(x) <= 0.2) {
+      printf("Absolute difference too large !\n");
+      return 1;
+    }
+    if (abs > 1.2e-7 && fabsf(x) <= 0.1) {
+      printf("Absolute difference too large !\n");
+      return 1;
+    }
+
+    if (rel > 1e-6 && fabsf(x) <= 0.2) {
+      printf("Relative difference too large !\n");
+      return 1;
+    }
+    if (rel > 4e-8 && fabsf(x) <= 0.1) {
+      printf("Relative difference too large !\n");
+      return 1;
+    }
+
+  }
+
+  printf("\nAll values are consistent\n");
+
+  return 0;
+}