I am trying to find the longest Euclidean distance to the origin for each point. Their positions are stored in two float
arrays, denoting each point's position on the Cartesian plane. I tried to speed up the function using SIMD (WITHOUT SVML, so no _mm256_hypot_ps
), however it still seemed a little bit slow. The code loop through the array by a chunk size of 8, on each iteration updates the maximum squared Euclidean distance and finally returns the sqrt of the maximum element.
#include <immintrin.h>
#include <math.h>
#include <stdio.h>
float dist(const float* x, const float* y, unsigned int len) {
if (len < 16) {
float d = 0.0f, a;
unsigned i = 0;
for (i = 0; i < len; i++) {
a = x[i] * x[i] + y[i] * y[i];
if (a > d) d = a;
}
return sqrt(d);
}
__m256 a = _mm256_loadu_ps(x), b = _mm256_loadu_ps(y);
a = _mm256_mul_ps(a, a);
b = _mm256_mul_ps(b, b);
__m256 max = _mm256_add_ps(a, b); // x[i]*x[i] + y[i]*y[i]
unsigned i;
for (i = 8; i + 8 < len; i += 8) {
a = _mm256_loadu_ps(x + i);
a = _mm256_mul_ps(a, a);
b = _mm256_loadu_ps(y + i);
b = _mm256_mul_ps(b, b);
max = _mm256_max_ps(
max, _mm256_add_ps(a, b)); // max(max[i], x[i]*x[i] + y[i]*y[i])
}
float d = 0.0f, t;
for (; i < len; ++i) {
t = x[i] * x[i] + y[i] * y[i];
if (t > d) d = t;
}
const __m128 hiQuad = _mm256_extractf128_ps(max, 1); // high 128 bits in max
const __m128 loQuad = _mm256_castps256_ps128(max); // low 128 bits in max
const __m128 maxQuad = _mm_max_ps(loQuad, hiQuad);
const __m128 hiDual = _mm_movehl_ps(maxQuad, maxQuad);
const __m128 loDual = maxQuad;
const __m128 maxDual = _mm_max_ps(loDual, hiDual);
const __m128 lo = maxDual;
const __m128 hi = _mm_shuffle_ps(maxDual, maxDual, 0x1);
const __m128 res = _mm_max_ss(lo, hi);
t = _mm_cvtss_f32(res);
if (t > d) d = t;
return sqrt(d);
}
The function will be called many times. I would like to speed it up further.