test igraph umap layout with avx2

umap.avx2.patch

 src/layout/umap.c | 58 +++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 58 insertions(+)

diff --git a/src/layout/umap.c b/src/layout/umap.c
index d1f334313..a9e40d7d7 100644
--- a/src/layout/umap.c
+++ b/src/layout/umap.c
@@ -28,6 +28,7 @@

 #include "layout/layout_internal.h"
 #include "core/interruption.h"
+#include <immintrin.h>

 #include <math.h>

@@ -720,6 +721,7 @@ static igraph_real_t igraph_i_umap_clip_force(igraph_real_t force, igraph_real_t
     return force > limit ? limit : (force < -limit ? -limit : force);
 }

+/*
 static igraph_real_t igraph_i_umap_attract(
         igraph_real_t dsq,
         igraph_real_t a,
@@ -727,7 +729,42 @@ static igraph_real_t igraph_i_umap_attract(
 {
     return - (2 * a * b * pow(dsq, b - 1.)) / (1. + a * pow(dsq, b));
 }
+*/
+
+
+#ifndef _mm256_extract_pd
+#define _mm256_extract_pd(R, I) (*((double*)(&R)+I))
+#endif
+
+__m256d _mm256_pow_pd(__m256d base, __m256d exponent) {
+    __m256d result = _mm256_set1_pd(1.0);
+    __m256d base_n = base;
+    __m256d exponent_n = exponent;
+    for (int i = 0; i < 4; i++) {
+        if (_mm256_extract_pd(exponent_n, i) >= 1)
+            result = _mm256_mul_pd(result, base_n);
+        base_n = _mm256_mul_pd(base_n, base);
+        exponent_n = _mm256_sub_pd(exponent_n, _mm256_set1_pd(1));
+    }
+    return result;
+}
+
+static igraph_real_t igraph_i_umap_attract(
+          igraph_real_t dsq,
+          igraph_real_t a,
+          igraph_real_t b)
+{
+    __m256d dsq_vec = _mm256_set1_pd(dsq);
+    __m256d a_vec = _mm256_set1_pd(a);
+    __m256d b_vec = _mm256_set1_pd(b);
+    __m256d temp1 = _mm256_mul_pd(a_vec, _mm256_pow_pd(dsq_vec, b_vec));
+    __m256d temp2 = _mm256_add_pd(temp1, _mm256_set1_pd(1));
+    __m256d temp3 = _mm256_mul_pd(_mm256_mul_pd(b_vec, _mm256_set1_pd(2)), _mm256_pow_pd(dsq_vec, _mm256_sub_pd(b_vec, _mm256_set1_pd(1))));
+    __m256d result = _mm256_div_pd(_mm256_mul_pd(temp3, a_vec), temp2);
+    return -_mm256_extract_pd(result, 0);
+}

+/*
 static igraph_real_t igraph_i_umap_repel(
         igraph_real_t dsq,
         igraph_real_t a,
@@ -737,6 +774,27 @@ static igraph_real_t igraph_i_umap_repel(

     return (2 * b) / (dsq_min + dsq) / (1. + a * pow(dsq, b));
 }
+*/
+
+
+
+static igraph_real_t igraph_i_umap_repel(igraph_real_t dsq, igraph_real_t a, igraph_real_t b) {
+    igraph_real_t dsq_min = UMAP_CORRECT_DISTANCE_REPULSION * UMAP_CORRECT_DISTANCE_REPULSION;
+
+    __m256d dsq_min_vec = _mm256_set1_pd(dsq_min);
+    __m256d dsq_vec = _mm256_set1_pd(dsq);
+    __m256d a_vec = _mm256_set1_pd(a);
+    __m256d b_vec = _mm256_set1_pd(b);
+    __m256d dsq_pow_b = _mm256_pow_pd(dsq_vec, b_vec);
+    __m256d dsq_pow_b_a = _mm256_mul_pd(dsq_pow_b, a_vec);
+    __m256d one_plus_dsq_pow_b_a = _mm256_add_pd(dsq_pow_b_a, _mm256_set1_pd(1.0));
+    __m256d dsq_plus_dsq_min = _mm256_add_pd(dsq_vec, dsq_min_vec);
+    __m256d two_b_vec = _mm256_set1_pd(2.0 * b);
+    __m256d result = _mm256_div_pd(two_b_vec, _mm256_div_pd(dsq_plus_dsq_min, one_plus_dsq_pow_b_a));
+    return _mm256_extract_pd(result,0);
+}
+
+

 static igraph_error_t igraph_i_umap_apply_forces(
         const igraph_t *graph,