diff --git a/math_intrinsics.h b/math_intrinsics.h
index 53eb3e8..3efb5be 100644
--- a/math_intrinsics.h
+++ b/math_intrinsics.h
@@ -52,6 +52,9 @@ extern "C" {
     // max error : 4.768371582e-07
     float32x4_t vcbrtq_f32(float32x4_t a);
 
+    // max error : 1.484901873e-07
+    float32x4_t vpowq_f32(float32x4_t x, float32x4_t y);
+
     #define __MATH__INTRINSICS__NEON__
 
 #else
@@ -93,6 +96,9 @@ extern "C" {
     // max error : 4.768371582e-07
     __m256 mm256_cbrt_ps(__m256 a);
 
+    // max error : 1.484901873e-07
+    __m256 mm256_pow_ps(__m256 x, __m256 y);
+
     #define __MATH__INTRINSICS__AVX__
 
 #endif
@@ -167,11 +173,22 @@ extern "C" {
     static inline simd_vectori simd_andnot_i(simd_vectori a, simd_vectori b) {return vbicq_s32(a, b);}
     static inline simd_vectori simd_cmp_eq_i(simd_vectori a, simd_vectori b) {return vceqq_s32(a, b);}
     static inline simd_vectori simd_cmp_gt_i(simd_vectori a, simd_vectori b) {return vcgtq_s32(a, b);}
+    static inline simd_vectori simd_min_i(simd_vectori a, simd_vectori b) {return vminq_s32(a, b);}
+    static inline simd_vectori simd_max_i(simd_vectori a, simd_vectori b) {return vmaxq_s32(a, b);}
+    static inline simd_vector simd_gather(const float* array, simd_vectori indices)
+    {
+        float tmp[4] = {array[indices[0]], array[indices[1]], array[indices[2]], array[indices[3]]};
+        return vld1q_f32(tmp);
+    }
 
     #define simd_asin vasinq_f32
     #define simd_atan vatanq_f32
     #define simd_sincos vsincosq_f32
     #define simd_sin vsinq_f32
+    #define simd_log vlogq_f32
+    #define simd_exp vexpq_f32
+    #define simd_log2 vlog2q_f32
+    #define simd_exp2 vexp2q_f32
 
 #else
     typedef __m256 simd_vector;
@@ -231,17 +248,25 @@ extern "C" {
     static inline simd_vectori simd_splat_i(int i) {return _mm256_set1_epi32(i);}
     static inline simd_vectori simd_splat_zero_i(void) {return _mm256_setzero_si256();}
     static inline simd_vectori simd_shift_left_i(simd_vectori a, int i) {return _mm256_slli_epi32(a, i);}
-    static inline simd_vectori simd_shift_right_i(simd_vectori a, int i) {return _mm256_srli_epi32(a, i);}
+    static inline simd_vectori simd_shift_right_i(simd_vectori a, int i) {return _mm256_srai_epi32(a, i);}
     static inline simd_vectori simd_and_i(simd_vectori a, simd_vectori b) {return _mm256_and_si256(a, b);}
     static inline simd_vectori simd_or_i(simd_vectori a, simd_vectori b) {return _mm256_or_si256(a, b);}
     static inline simd_vectori simd_andnot_i(simd_vectori a, simd_vectori b) {return _mm256_andnot_si256(b, a);}
     static inline simd_vectori simd_cmp_eq_i(simd_vectori a, simd_vectori b) {return _mm256_cmpeq_epi32(a, b);}
     static inline simd_vectori simd_cmp_gt_i(simd_vectori a, simd_vectori b) {return _mm256_cmpgt_epi32(a, b);}
+    static inline simd_vectori simd_min_i(simd_vectori a, simd_vectori b) {return _mm256_min_epi32(a, b);}
+    static inline simd_vectori simd_max_i(simd_vectori a, simd_vectori b) {return _mm256_max_epi32(a, b);}
+    static inline simd_vector simd_gather(const float* array, simd_vectori indices) {return _mm256_i32gather_ps(array, indices, 4);}
+
 
     #define simd_asin mm256_asin_ps
     #define simd_atan mm256_atan_ps
     #define simd_sincos mm256_sincos_ps
     #define simd_sin mm256_sin_ps
+    #define simd_exp mm256_exp_ps
+    #define simd_log mm256_log_ps
+    #define simd_exp2 mm256_exp2_ps
+    #define simd_log2 mm256_log2_ps
 
 #endif
 
@@ -309,6 +334,10 @@ static inline simd_vector simd_sign(simd_vector a)
     return simd_select(result, simd_splat( 1.f), simd_cmp_gt(a, simd_splat_zero()));
 }
 
+static inline simd_vectori simd_select_i(simd_vectori a, simd_vectori b, simd_vectori mask) { return simd_or_i(simd_andnot_i(a, mask), simd_and_i(b, mask));}
+static inline simd_vectori simd_neg_i(simd_vectori a){return simd_sub_i(simd_splat_zero_i(), a);}
+
+
 //----------------------------------------------------------------------------------------------------------------------
 // based on http://gruntthepeon.free.fr/ssemath/
 #ifdef __MATH__INTRINSICS__NEON__
@@ -810,6 +839,195 @@ __m256 mm256_cos_ps(__m256 x)
     return x;
 }
 
+static inline simd_vector reduc(simd_vector x) {return simd_mul(simd_splat(0.0625f),  simd_floor( simd_mul(simd_splat(16.f),x)));}
+
+//----------------------------------------------------------------------------------------------------------------------
+// based on https://github.com/jeremybarnes/cephes/blob/master/single/powf.c
+#ifdef __MATH__INTRINSICS__NEON__
+    float32x4_t vpowq_f32(float32x4_t x, float32x4_t y)
+#else
+    __m256 mm256_pow_ps(__m256 x, __m256 y)
+#endif
+{
+    simd_vector x_equals_zero = simd_cmp_eq(x, simd_splat_zero());
+    simd_vector y_equals_zero = simd_cmp_eq(y, simd_splat_zero());
+    simd_vector non_integer_power = simd_cmp_neq(y, simd_floor(y));
+    simd_vector return_zero = simd_andnot(x_equals_zero, y_equals_zero);
+    simd_vector return_one = simd_and(x_equals_zero, y_equals_zero);
+    simd_vector return_nan = simd_and(simd_cmp_lt(x, simd_splat_zero()), non_integer_power);
+
+#ifdef __MATH_INTRINSINCS_FAST__
+    simd_vector z = simd_exp2(simd_mul(y, simd_log2(x)));
+#else
+    // 2^(-i/16) The decimal values are rounded to 24-bit precision
+    static float A[] =
+    {
+        1.00000000000000000000E0f,
+        9.57603275775909423828125E-1f,
+        9.17004048824310302734375E-1f,
+        8.78126084804534912109375E-1f,
+        8.40896427631378173828125E-1f,
+        8.05245161056518554687500E-1f,
+        7.71105408668518066406250E-1f,
+        7.38413095474243164062500E-1f,
+        7.07106769084930419921875E-1f,
+        6.77127778530120849609375E-1f,
+        6.48419797420501708984375E-1f,
+        6.20928883552551269531250E-1f,
+        5.94603538513183593750000E-1f,
+        5.69394290447235107421875E-1f,
+        5.45253872871398925781250E-1f,
+        5.22136867046356201171875E-1f,
+        5.00000000000000000000E-1f
+    };
+
+    // continuation, for even i only 2^(i/16)  =  A[i] + B[i/2]
+    static float B[] =
+    {
+        0.00000000000000000000E0f,
+        -5.61963907099083340520586E-9f,
+        -1.23776636307969995237668E-8f,
+        4.03545234539989593104537E-9f,
+        1.21016171044789693621048E-8f,
+        -2.00949968760174979411038E-8f,
+        1.89881769396087499852802E-8f,
+        -6.53877009617774467211965E-9f,
+        0.00000000000000000000E0f
+    };
+
+    // 1 / A[i] The decimal values are full precision
+    static float Ainv[] =
+    {
+        1.00000000000000000000000E0f,
+        1.04427378242741384032197E0f,
+        1.09050773266525765920701E0f,
+        1.13878863475669165370383E0f,
+        1.18920711500272106671750E0f,
+        1.24185781207348404859368E0f,
+        1.29683955465100966593375E0f,
+        1.35425554693689272829801E0f,
+        1.41421356237309504880169E0f,
+        1.47682614593949931138691E0f,
+        1.54221082540794082361229E0f,
+        1.61049033194925430817952E0f,
+        1.68179283050742908606225E0f,
+        1.75625216037329948311216E0f,
+        1.83400808640934246348708E0f,
+        1.91520656139714729387261E0f,
+        2.00000000000000000000000E0f
+    };
+
+    simd_vector neg_x = simd_andnot(simd_cmp_lt(x, simd_splat_zero()), non_integer_power);
+
+    x = simd_select(x, simd_neg(x), neg_x);
+
+    // separate significand from exponent
+    simd_vector e;
+    x = simd_frexp(x, &e);
+
+    // find significand in antilog table A[]
+    simd_vectori i = simd_splat_i(1);
+    i = simd_select_i(i, simd_splat_i(9), simd_cast_from_float(simd_cmp_le(x, simd_splat(A[9]))));
+    simd_vectori i_plus_4 = simd_add_i(i, simd_splat_i(4));
+    i = simd_select_i(i, i_plus_4, simd_cast_from_float(simd_cmp_le(x, simd_gather(A, i_plus_4))));
+    simd_vectori i_plus_2 = simd_add_i(i, simd_splat_i(2));
+    i = simd_select_i(i, i_plus_2, simd_cast_from_float(simd_cmp_le(x, simd_gather(A, i_plus_2))));
+    i = simd_select_i(i, simd_splat_i(-1), simd_cast_from_float(simd_cmp_ge(x, simd_splat(A[1]))));
+    i = simd_add_i(i, simd_splat_i(1));
+
+    // Find (x - A[i])/A[i]
+    // in order to compute log(x/A[i]):
+    // log(x) = log( a x/a ) = log(a) + log(x/a)
+    // log(x/a) = log(1+v),  v = x/a - 1 = (x-a)/a
+    x = simd_sub(x, simd_gather(A, i));
+    x = simd_sub(x, simd_gather(B, simd_shift_right_i(i, 1)));
+    x = simd_mul(x, simd_gather(Ainv, i));
+
+    // rational approximation for log(1+v):
+    // log(1+v)  =  v  -  0.5 v^2  +  v^3 P(v)
+    // Theoretical relative error of the approximation is 3.5e-11
+    // on the interval 2^(1/16) - 1  > v > 2^(-1/16) - 1
+    simd_vector z = simd_mul(x, x);
+    simd_vector w = simd_polynomial4(x, (float[]){-0.1663883081054895f, 0.2003770364206271f, -0.2500006373383951f, 0.3333331095506474f});
+    w = simd_mul(w, x);
+    w = simd_fmad(w, z, simd_mul(simd_splat(-.5f), z));
+
+    // Convert to base 2 logarithm: multiply by log2(e)
+    simd_vector LOG2EA = simd_splat(0.44269504088896340736F);
+    w = simd_fmad(w, LOG2EA, w);
+
+    // Note x was not yet added in to above rational approximation,
+    // so do it now, while multiplying by log2(e).
+    z = simd_fmad(x, LOG2EA, w);
+    z = simd_add(z, x);
+
+    // Compute exponent term of the base 2 logarithm.
+    w = simd_neg(simd_convert_from_int(i));
+    w = simd_fmad(w, simd_splat(0.0625f), e);
+
+    // Multiply base 2 log by y, in extended precision. 
+    // separate y into large part ya and small part yb less than 1/16
+    simd_vector ya = reduc(y);
+    simd_vector yb = simd_sub(y, ya);
+
+    simd_vector W = simd_fmad(z, y, simd_mul(w, yb));
+    simd_vector Wa = reduc(W);
+    simd_vector Wb = simd_sub(W, Wa);
+
+    W = simd_fmad(w, ya, Wa);
+    Wa = reduc(W);
+    simd_vector u = simd_sub(W, Wa);
+
+    W = simd_add(Wb, u);
+    Wb = reduc(W);
+    w = simd_mul(simd_splat(16.f), simd_add(Wa, Wb));
+
+    return_zero = simd_or(return_zero, simd_cmp_lt(w, simd_splat(-2400.0f)));
+
+    e = w;
+    Wb = simd_sub(W, Wb);
+
+    simd_vector gt_zero = simd_cmp_gt(Wb, simd_splat_zero());
+    e = simd_select(e, simd_add(e, simd_splat(1.f)), gt_zero);
+    Wb = simd_select(Wb, simd_sub(Wb, simd_splat(0.0625f)), gt_zero);
+
+    // Now the product y * log2(x)  =  Wb + e/16.0.
+    // Compute base 2 exponential of Wb,
+    // where -0.0625 <= Wb <= 0.
+    // Theoretical relative error of the approximation is 2.8e-12.
+    //  z  =  2**Wb - 1
+    z = simd_polynomial4(Wb, (float[]) {9.416993633606397E-003f, 5.549356188719141E-002f, 2.402262883964191E-001f, 6.931471791490764E-001f});
+    z = simd_mul(z, Wb);
+
+    simd_vector neg_e = simd_cmp_lt(e, simd_splat_zero());
+    simd_vectori int_e = simd_convert_from_float(e);
+
+    simd_vectori i0 = simd_neg_i(simd_shift_right_i(simd_neg_i(int_e), 4));
+    simd_vectori i1 = simd_add_i(simd_shift_right_i(int_e, 4), simd_splat_i(1));
+    i = simd_select_i(i1, i0, simd_cast_from_float(neg_e));
+
+    int_e = simd_sub_i(simd_shift_left_i(i, 4), int_e);
+
+    // clamp int_e to avoid reading data out of the array 
+    int_e = simd_min_i(int_e, simd_splat_i(16));
+    int_e = simd_max_i(int_e, simd_splat_i(0));
+
+    w = simd_gather(A, int_e);
+    z = simd_fmad(w, z, w);     // 2^-e * ( 1 + (2^Hb-1) )
+    z = simd_ldexp(z, simd_convert_from_int(i));
+
+    // For negative x, find out if the integer exponent is odd or even.
+    w = simd_mul(simd_splat(2.f), simd_floor(simd_mul(simd_splat(.5f), w)));
+    z = simd_select(z, simd_neg(z), simd_and(neg_x, simd_cmp_neq(w, y)));
+#endif
+
+    z = simd_andnot(z, return_zero);
+    z = simd_select(z, simd_splat(1.f), return_one);
+    z = simd_or(z, return_nan);
+
+    return z;
+}
+
 #endif  // __MATH__INTRINSICS__IMPLEMENTATION__
 
 
diff --git a/tests/benchmark.c b/tests/benchmark.c
index 809dad9..9af6a19 100644
--- a/tests/benchmark.c
+++ b/tests/benchmark.c
@@ -7,17 +7,23 @@
 
 #include "../math_intrinsics.h"
 
+//----------------------------------------------------------------------------------------------------------------------
+// functions pointer definition
 typedef float (*reference_function)(float);
+typedef float (*reference_function2)(float, float);
 #ifdef __MATH__INTRINSICS__AVX__
     typedef __m256 (*approximation_function)(__m256);
+    typedef __m256 (*approximation_function2)(__m256, __m256);
     #define simd_vector_width (8)
 #else
     typedef float32x4_t (*approximation_function)(float32x4_t);
+    typedef float32x4_t (*approximation_function2)(float32x4_t, float32x4_t);
     #define simd_vector_width (4)
 #endif
 
 #define NUM_ITERATIONS (200000000)
 
+//----------------------------------------------------------------------------------------------------------------------
 int benchmark(approximation_function function, reference_function reference, const char* name)
 {
     float init_array[simd_vector_width];
@@ -59,9 +65,9 @@ int benchmark(approximation_function function, reference_function reference, con
 
     float simd_time = stm_ms(stm_since(start));
 
-    printf(".%s:\t %5.2fms", name, simd_time);
+    printf(".%s:\t %05.2fms", name, simd_time);
 
-    float total = simd_time;
+    float total = simd_time / 1000.f;
 
     start = stm_now();
 
@@ -72,12 +78,83 @@ int benchmark(approximation_function function, reference_function reference, con
 
     float clib_time = stm_ms(stm_since(start));
     
-    printf("\tc std func: %5.2fms\tratio: %2.2fx\n", clib_time, clib_time/simd_time);
+    printf("\tc std func: %05.2fms\tratio: %2.2fx\n", clib_time, clib_time/simd_time);
 
     return output;
 }
 
+//----------------------------------------------------------------------------------------------------------------------
+int benchmark2(approximation_function2 function, reference_function2 reference, const char* name)
+{
+    float array_x[simd_vector_width];
+    float array_y[simd_vector_width];
+    uint64_t start = 0;
+    int output = 0;
+
+    for(uint32_t i=0; i<simd_vector_width; ++i)
+    {
+        array_x[i] = (float) (i) / (float) (simd_vector_width);
+        array_y[i] = (float) (i); 
+    }
+
+#ifdef __MATH__INTRINSICS__AVX__
+    __m256 step = _mm256_set1_ps(FLT_EPSILON);
+    __m256 v_x = _mm256_loadu_ps(array_x);
+    __m256 v_y = _mm256_loadu_ps(array_y);
+    __m256 result = _mm256_setzero_ps();
+
+     start = stm_now();
+
+    for(uint32_t i=0; i<(NUM_ITERATIONS / simd_vector_width); ++i)
+    {
+        result = _mm256_add_ps(result, function(v_x, v_y));
+        v_x = _mm256_add_ps(v_x, step);
+    }
+
+    output = _mm256_cvtss_f32(result);
+#else
+    float32x4_t step = vdupq_n_f32(FLT_EPSILON);
+    float32x4_t v_x = vld1q_f32(array_x);
+    float32x4_t v_y = vld1q_f32(array_y);
+    float32x4_t result = vdupq_n_f32(0.f);
+
+     start = stm_now();
+
+    for(uint32_t i=0; i<(NUM_ITERATIONS / simd_vector_width); ++i)
+    {
+        result = vaddq_f32(result, function(v_x, v_y));
+        v_x = vaddq_f32(v_x, step);
+    }
+
+    output = vgetq_lane_f32(result, 0);
+#endif
+
+    float simd_time = stm_ms(stm_since(start));
+
+    printf(".%s:\t %05.2fms", name, simd_time);
+
+    float total = simd_time / 1000.f;
+
+    start = stm_now();
+
+    float x = 0.f;
+
+    for(uint32_t i=0; i<NUM_ITERATIONS; ++i)
+    {
+        total += reference(x, x*2.f);
+        x += 0.001f;
+    }
+
+    output += total;
+
+    float clib_time = stm_ms(stm_since(start));
+    
+    printf("\tc std func: %05.2fms\tratio: %2.2fx\n", clib_time, clib_time/simd_time);
+
+    return output;
+}
 
+//----------------------------------------------------------------------------------------------------------------------
 int main(int argc, char * argv[])
 {
     stm_setup();
@@ -94,6 +171,7 @@ int main(int argc, char * argv[])
     output += benchmark(mm256_acos_ps, acosf, "mm256_acos_ps");
     output += benchmark(mm256_asin_ps, asinf, "mm256_asin_ps");
     output += benchmark(mm256_atan_ps, atanf, "mm256_atan_ps");
+    output += benchmark2(mm256_atan2_ps, atan2f, "mm256_atan2_ps");
     output += benchmark(mm256_cbrt_ps, cbrtf, "mm256_cbrt_ps");
     output += benchmark(mm256_cos_ps, cosf, "mm256_cos_ps");
     output += benchmark(mm256_sin_ps, sinf, "mm256_sin_ps");
@@ -101,10 +179,12 @@ int main(int argc, char * argv[])
     output += benchmark(mm256_exp2_ps, exp2f, "mm256_exp2_ps");
     output += benchmark(mm256_log_ps, logf, "mm256_log_ps");
     output += benchmark(mm256_log2_ps, log2f, "mm256_log2_ps");
+    output += benchmark2(mm256_pow_ps, powf, "mm256_pow_ps");
 #else
     output += benchmark(vacosq_f32, acosf, "vacosq_f32");
     output += benchmark(vasinq_f32, asinf, "vasinq_f32");
     output += benchmark(vatanq_f32, atanf, "vatanq_f32");
+    output += benchmark2(vatan2q_f32, atan2f, "vatan2q_f32");
     output += benchmark(vcbrtq_f32, cbrtf, "vcbrtq_f32");
     output += benchmark(vcosq_f32, cosf, "vcosq_f32");
     output += benchmark(vsinq_f32, sinf, "vsinq_f32");
@@ -112,6 +192,7 @@ int main(int argc, char * argv[])
     output += benchmark(vexp2q_f32, exp2f, "vexp2q_f32");
     output += benchmark(vlogq_f32, logf, "vlogq_f32");
     output += benchmark(vlog2q_f32, log2f, "vlog2q_f32");
+    output += benchmark2(vpowq_f32, powf, "vpowq_f32");
 #endif
 
     printf("\n%d\n", output);
diff --git a/tests/test.c b/tests/test.c
index 9c02ff0..674fe43 100644
--- a/tests/test.c
+++ b/tests/test.c
@@ -17,11 +17,14 @@
 //----------------------------------------------------------------------------------------------------------------------
 // functions pointer definition
 typedef float (*reference_function)(float);
+typedef float (*reference_function2)(float, float);
 #ifdef __MATH__INTRINSICS__AVX__
     typedef __m256 (*approximation_function)(__m256);
+    typedef __m256 (*approximation_function2)(__m256, __m256);
     #define simd_vector_width (8)
 #else
     typedef float32x4_t (*approximation_function)(float32x4_t);
+    typedef float32x4_t (*approximation_function2)(float32x4_t, float32x4_t);
     #define simd_vector_width (4)
 #endif
 
@@ -90,6 +93,76 @@ TEST generic_test(reference_function ref, approximation_function approx, float r
     PASS();
 }
 
+//----------------------------------------------------------------------------------------------------------------------
+// generic unit test with 2 arguments
+TEST generic_test2(reference_function2 ref, approximation_function2 approx, float range_min, float range_max, float epsilon, uint32_t num_elements, bool relative_error, const char* name)
+{
+    float* x = (float*) malloc(num_elements * sizeof(float));
+    float* y = (float*) malloc(num_elements * sizeof(float));
+    float* result = (float*) malloc(num_elements * sizeof(float));
+    float step = ((range_max - range_min) / (float) (num_elements-1));
+    uint32_t num_vectors = num_elements / simd_vector_width;
+
+    for(uint32_t i=0; i<num_elements; ++i)
+    {
+        x[i] = (step * (float)(i)) + range_min;
+        y[i] = sinf(x[i]);
+        result[i] = ref(x[i], y[i]);
+    }
+
+#ifdef __MATH__INTRINSICS__AVX__
+    __m256 v_epsilon = _mm256_set1_ps(epsilon);
+    __m256 v_max_error = _mm256_setzero_ps();
+
+    for(uint32_t i=0; i<num_vectors; ++i)
+    {
+        __m256 v_x = _mm256_loadu_ps(x + i * simd_vector_width);
+        __m256 v_y = _mm256_loadu_ps(y + i * simd_vector_width);
+        __m256 v_result = _mm256_loadu_ps(result+ i * simd_vector_width);
+        __m256 v_approx = approx(v_x, v_y);
+        __m256 v_error = _mm256_and_ps(_mm256_sub_ps(v_approx, v_result), _mm256_castsi256_ps(_mm256_set1_epi32(0x7FFFFFFF)));
+
+        if (relative_error)
+            v_error = _mm256_div_ps(v_error, v_result);
+        
+        ASSERT(_mm256_movemask_ps(_mm256_cmp_ps(v_error, v_epsilon, _CMP_LE_OQ)) == 0xff);
+        v_max_error = _mm256_max_ps(v_max_error, v_error);
+    }
+
+    v_max_error = _mm256_max_ps(v_max_error, _mm256_permute_ps(v_max_error, _MM_SHUFFLE(2, 1, 0, 3)));
+    v_max_error = _mm256_max_ps(v_max_error, _mm256_permute_ps(v_max_error, _MM_SHUFFLE(1, 0, 3, 2)));
+    v_max_error = _mm256_max_ps(v_max_error, _mm256_permute2f128_ps(v_max_error, v_max_error, 1));
+
+    printf("%s\t max error : %.*e\n", name, FLT_DECIMAL_DIG, _mm256_cvtss_f32(v_max_error));
+#else
+    float32x4_t v_epsilon = vdupq_n_f32(epsilon);
+    float32x4_t v_max_error = vdupq_n_f32(0.f);
+
+    for(uint32_t i=0; i<num_vectors; ++i)
+    {
+        float32x4_t v_x = vld1q_f32(x + i * simd_vector_width);
+        float32x4_t v_y = vld1q_f32(y + i * simd_vector_width);
+        float32x4_t v_result = vld1q_f32(result+ i * simd_vector_width);
+        float32x4_t v_approx = approx(v_x, v_y);
+        float32x4_t v_error = vabsq_f32(vsubq_f32(v_approx, v_result));
+
+        if (relative_error)
+            v_error = vdivq_f32(v_error, v_result);
+        
+        ASSERT(vminvq_u32(vcleq_f32(v_error, v_epsilon)) == UINT32_MAX);
+        v_max_error = vmaxq_f32(v_max_error, v_error);
+    }
+
+    printf("%s\t max error : %.*e\n", name, FLT_DECIMAL_DIG, vmaxvq_f32(v_max_error));
+#endif
+
+    free(x);
+    free(y);
+    free(result);
+    
+    PASS();
+}
+
 //----------------------------------------------------------------------------------------------------------------------
 TEST value_expected(float input, float target, approximation_function function)
 {
@@ -122,38 +195,23 @@ TEST nan_expected(float input, approximation_function function)
     PASS();
 }
 
-
-//----------------------------------------------------------------------------------------------------------------------
-float atan2_angle(float angle) {return atan2f(sinf(angle) * (angle + 1.f), cosf(angle) * (angle + 1.f));}
-
-#ifdef __MATH__INTRINSICS__AVX__
-__m256 simd_atan2(__m256 angle) 
-{
-    __m256 angle_plus_one = _mm256_add_ps(angle, _mm256_set1_ps(1.f));
-    return mm256_atan2_ps(_mm256_mul_ps(mm256_cos_ps(angle), angle_plus_one), _mm256_mul_ps(mm256_sin_ps(angle), angle_plus_one));}
-#else
-float32x4_t simd_atan2(float32x4_t angle) 
-{
-    float32x4_t angle_plus_one = vaddq_f32(angle, vdupq_n_f32(1.f));
-    return vatan2q_f32(vmulq_f32(vcosq_f32(angle), angle_plus_one), vmulq_f32(vsinq_f32(angle), angle_plus_one));
-}
-#endif
-
-
 #define NUM_SAMPLES (1024)
 
+static const float cbrt_threshold = 2.e-07f;
+
 #ifdef __MATH_INTRINSINCS_FAST__
 static const float trigo_threshold = 1.e-06f;
 static const float exp_threshold = 3.e-07f;
-static const float arc_theshold = 1.e-04f;
-static const float cbrt_threshold = 2.e-04f;
+static const float arc_threshold = 1.e-04f;
+static const float pow_threshold = 1.e-05f;
 #else
 static const float trigo_threshold = FLT_EPSILON;
 static const float exp_threshold = 2.e-07f;
-static const float arc_theshold = 1.e-06f;
-static const float cbrt_threshold = 2.e-07f;
+static const float arc_threshold = 1.e-06f;
+static const float pow_threshold = 2.e-07f;
 #endif
 
+float atan2_xy(float x, float y) {return atan2f(y, x);}
 
 SUITE(trigonometry)
 {
@@ -162,21 +220,21 @@ SUITE(trigonometry)
 #ifdef __MATH__INTRINSICS__AVX__
     RUN_TESTp(generic_test, sinf, mm256_sin_ps, -10.f, 10.f, trigo_threshold, NUM_SAMPLES, false, "mm256_sin_ps");
     RUN_TESTp(generic_test, cosf, mm256_cos_ps, -10.f, 10.f, trigo_threshold, NUM_SAMPLES, false, "mm256_cos_ps");
-    RUN_TESTp(generic_test, acosf, mm256_acos_ps, -1.f, 1.f, arc_theshold, NUM_SAMPLES, false, "mm256_acos_ps");
-    RUN_TESTp(generic_test, asinf, mm256_asin_ps, -1.f, 1.f, arc_theshold, NUM_SAMPLES, false, "mm256_asin_ps");
-    RUN_TESTp(generic_test, atanf, mm256_atan_ps, -10.f, 10.f, arc_theshold, NUM_SAMPLES, false, "mm256_atan_ps");
+    RUN_TESTp(generic_test, acosf, mm256_acos_ps, -1.f, 1.f, arc_threshold, NUM_SAMPLES, false, "mm256_acos_ps");
+    RUN_TESTp(generic_test, asinf, mm256_asin_ps, -1.f, 1.f, arc_threshold, NUM_SAMPLES, false, "mm256_asin_ps");
+    RUN_TESTp(generic_test, atanf, mm256_atan_ps, -10.f, 10.f, arc_threshold, NUM_SAMPLES, false, "mm256_atan_ps");
 
     // this task fails on linux and I don't have this OS to debug
     #if !defined(__linux__)
-        RUN_TESTp(generic_test, atan2_angle, simd_atan2, 0.f, 6.28318530f, 100.0f, 32768, false, "mm256_atan2_ps");
+        RUN_TESTp(generic_test2, atan2_xy, mm256_atan2_ps, 0.f, 1000.f, 3.e-07f, NUM_SAMPLES, false, "mm256_atan2_ps");
     #endif
 #else
     RUN_TESTp(generic_test, sinf, vsinq_f32, -10.f, 10.f, trigo_threshold, NUM_SAMPLES, false, "vsinq_f32");
     RUN_TESTp(generic_test, cosf, vcosq_f32, -10.f, 10.f, trigo_threshold, NUM_SAMPLES, false, "vcosq_f32");
-    RUN_TESTp(generic_test, acosf, vacosq_f32, -1.f, 1.f, arc_theshold, NUM_SAMPLES, false, "vacosq_f32");
-    RUN_TESTp(generic_test, asinf, vasinq_f32, -1.f, 1.f, arc_theshold, NUM_SAMPLES, false, "vasinq_f32");
-    RUN_TESTp(generic_test, atanf, vatanq_f32, -10.f, 10.f, arc_theshold, NUM_SAMPLES, false, "vatanq_f32");
-    RUN_TESTp(generic_test, atan2_angle, simd_atan2, 0.f, 6.28318530f, arc_theshold, 32768, false, "vatan2q_f32");
+    RUN_TESTp(generic_test, acosf, vacosq_f32, -1.f, 1.f, arc_threshold, NUM_SAMPLES, false, "vacosq_f32");
+    RUN_TESTp(generic_test, asinf, vasinq_f32, -1.f, 1.f, arc_threshold, NUM_SAMPLES, false, "vasinq_f32");
+    RUN_TESTp(generic_test, atanf, vatanq_f32, -10.f, 10.f, arc_threshold, NUM_SAMPLES, false, "vatanq_f32");
+    RUN_TESTp(generic_test2, atan2_xy, vatan2q_f32, 0.f, 1000.f, 3.e-07f, NUM_SAMPLES, false, "vatan2q_f32");
 #endif
 }
 
@@ -189,12 +247,14 @@ SUITE(exponentiation)
     RUN_TESTp(generic_test, expf, mm256_exp_ps, -87.f, 87.f, exp_threshold, NUM_SAMPLES, true, "mm256_exp_ps");
     RUN_TESTp(generic_test, exp2f, mm256_exp2_ps, -126.f, 126.f, exp_threshold, NUM_SAMPLES, true, "mm256_exp2");
     RUN_TESTp(generic_test, cbrtf, mm256_cbrt_ps, -1000.f, 1000.f, cbrt_threshold, 4096, true, "mm256_cbrt_ps");
+    RUN_TESTp(generic_test2, powf, mm256_pow_ps, 0.f, 100000.f, pow_threshold, 32768, true, "mm256_pow_ps");
 #else
     RUN_TESTp(generic_test, logf, vlogq_f32, FLT_EPSILON, 1.e20f, 1.e-07f, 32768, true, "vlogq_f32");
     RUN_TESTp(generic_test, log2f, vlog2q_f32, FLT_EPSILON, 1.e20f, 3.e-07f, 32768, true, "vlog2q_f32");
     RUN_TESTp(generic_test, expf, vexpq_f32, -87.f, 87.f, exp_threshold, NUM_SAMPLES, true, "vexpq_f32");
     RUN_TESTp(generic_test, exp2f, vexp2q_f32, -126.f, 126.f, exp_threshold, NUM_SAMPLES, true, "vexp2q_f32");
     RUN_TESTp(generic_test, cbrtf, vcbrtq_f32, -1000.f, 1000.f, cbrt_threshold, 4096, true, "vcbrtq_f32");
+    RUN_TESTp(generic_test2, powf, vpowq_f32, 0.f, 100000.f, pow_threshold, 32768, true, "vpowq_f32");
 #endif
 }