added simd_exp function and test

simd_math.c

@@ -1,7 +1,7 @@
 #include "simd_math.h"
 
 //-----------------------------------------------------------------------------
-// based on https://github.com/RJVB/sse_mathfun/blob/master/sse_mathfun.h
+// based on http://gruntthepeon.free.fr/ssemath/
 simd_vector simd_log(simd_vector x)
 {
     simd_vector one = simd_splat(1.f);
@@ -56,3 +56,53 @@ simd_vector simd_log(simd_vector x)
     return x;
 }
 
+//-----------------------------------------------------------------------------
+// based on http://gruntthepeon.free.fr/ssemath/
+simd_vector simd_exp(simd_vector x)
+{
+    simd_vector tmp = simd_splat_zero();
+    simd_vector fx;
+    simd_vector one = simd_splat(1.f);
+
+    x = simd_min(x, simd_splat(88.3762626647949f));
+    x = simd_max(x, simd_splat(-88.3762626647949f));
+
+    // express exp(x) as exp(g + n*log(2))
+    fx = simd_fmad(x, simd_splat(1.44269504088896341f), simd_splat(0.5f));
+    tmp = simd_floor(fx);
+
+    // if greater, substract 1
+    simd_vector mask = simd_cmp_gt(tmp, fx);
+    mask = simd_and(mask, one);
+    fx = simd_sub(tmp, mask);
+
+    tmp = simd_mul(fx, simd_splat(0.693359375f));
+    simd_vector z = simd_mul(fx, simd_splat(-2.12194440e-4f));
+    x = simd_sub(x, tmp);
+    x = simd_sub(x, z);
+    z = simd_mul(x, x);
+    simd_vector y = simd_splat(1.9875691500E-4f);
+    y = simd_fmad(y, x, simd_splat(1.3981999507E-3f));
+    y = simd_fmad(y, x, simd_splat(8.3334519073E-3f));
+    y = simd_fmad(y, x, simd_splat(4.1665795894E-2f));
+    y = simd_fmad(y, x, simd_splat(1.6666665459E-1f));
+    y = simd_fmad(y, x, simd_splat(5.0000001201E-1f));
+    y = simd_fmad(y, z, x);
+    y = simd_add(y, one);
+
+#if defined(SIMD_NEON_IMPLEMENTATION)
+    int32x4_t emm0 = vcvtq_s32_f32(fx);
+    emm0 = vaddq_s32(emm0, vdupq_n_s32(0x7f));
+    emm0 = vshlq_s32(emm0, vdupq_n_s32(23));
+    simd_vector pow2n = vreinterpretq_s32_f32(emm0);
+#elif defined(SIMD_AVX_IMPLEMENTATION)
+    __m256i emm0 = _mm256_cvtps_epi32(fx);
+    emm0 = _mm256_add_epi32(emm0, _mm256_set1_epi32(0x7f));
+    emm0 = _mm256_slli_epi32(emm0, 23);
+    simd_vector pow2n = _mm256_castsi256_ps(emm0);
+#endif
+
+    y = simd_mul(y, pow2n);
+    return y;
+}
+

simd_math.h

@@ -14,6 +14,7 @@ static simd_vector simd_acos(simd_vector x);
 static simd_vector simd_cos(simd_vector a);
 static simd_vector simd_atan2(simd_vector x, simd_vector y);
        simd_vector simd_log(simd_vector x);
+       simd_vector simd_exp(simd_vector x);
 static simd_vector simd_approx_cos(simd_vector a);
 static simd_vector simd_approx_sin(simd_vector a);
 static simd_vector simd_approx_exp(simd_vector x);

tests/test_simd.c

@@ -551,6 +551,7 @@ int main(int argc, char * argv[])
     RUN_TEST(sort);
     RUN_SUITE(horizontal);
     RUN_SUITE(trigonometry);
+    RUN_SUITE(exponentiation);
     RUN_SUITE(color_space);
     RUN_SUITE(arithmetic);
     RUN_SUITE(sqrt_and_rcp);

tests/test_simd_math.h

@@ -243,7 +243,38 @@ TEST logarithm(void)
         max_error = simd_max(max_error, relative_error);
     }
 
-    printf("simd_log max error : %f\n", simd_hmax(max_error));
+    printf(".simd_log max error : %f\n", simd_hmax(max_error));
+
+    PASS();
+}
+
+TEST exponential(void)
+{
+    float array[NUM_ELEMENTS];
+    float result[NUM_ELEMENTS];
+    float step = 174.f / (float) (NUM_ELEMENTS);
+
+    for(int i=0; i<NUM_ELEMENTS; ++i)
+    {
+        array[i] = (step * (float)i) - 87.f;
+        result[i] = expf(array[i]);
+    }
+
+    simd_vector epsilon = simd_splat(FLT_EPSILON);
+    simd_vector max_error = simd_splat_zero();
+
+    for(int i=0; i<NUM_VECTORS; ++i)
+    {
+        simd_vector v_array = simd_load_offset(array, i);
+        simd_vector v_result = simd_load_offset(result, i);
+        simd_vector v_approx = simd_exp(v_array);
+
+        simd_vector relative_error = simd_div(simd_abs_diff(v_approx, v_result), v_result);
+        ASSERT(simd_all(simd_cmp_lt(relative_error, epsilon)));
+        max_error = simd_max(max_error, relative_error);
+    }
+
+    printf("simd_exp max error : %f\n", simd_hmax(max_error));
 
     PASS();
 }
@@ -335,8 +366,13 @@ SUITE(trigonometry)
     RUN_TEST(approx_arcos);
     RUN_TEST(arctan);
     RUN_TEST(arctan2);
-    RUN_TEST(approx_exp);
+}
+
+SUITE(exponentiation)
+{
     RUN_TEST(logarithm);
+    RUN_TEST(exponential);
+    RUN_TEST(approx_exp);
 }
 
 SUITE(color_space)