As @Mysticial says in the comments above, do the compare and sum vertically and then just sum horizontally at the end of the main loop:
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <emmintrin.h>
// reference implementation
int fast_compare_ref(const char *s, const char *t, int length)
{
int result = 0;
int i;
for (i = 0; i < length; ++i)
{
if (s[i] == t[i])
result++;
}
return result;
}
// optimised implementation
int fast_compare(const char *s, const char *t, int length)
{
int result = 0;
int i;
__m128i vsum = _mm_set1_epi32(0);
for (i = 0; i < length - 15; i += 16)
{
__m128i vs, vt, v, vh, vl, vtemp;
vs = _mm_loadu_si128((__m128i *)&s[i]); // load 16 chars from input
vt = _mm_loadu_si128((__m128i *)&t[i]);
v = _mm_cmpeq_epi8(vs, vt); // compare
vh = _mm_unpackhi_epi8(v, v); // unpack compare result into 2 x 8 x 16 bit vectors
vl = _mm_unpacklo_epi8(v, v);
vtemp = _mm_madd_epi16(vh, vh); // accumulate 16 bit vectors into 4 x 32 bit partial sums
vsum = _mm_add_epi32(vsum, vtemp);
vtemp = _mm_madd_epi16(vl, vl);
vsum = _mm_add_epi32(vsum, vtemp);
}
// get sum of 4 x 32 bit partial sums
vsum = _mm_add_epi32(vsum, _mm_srli_si128(vsum, 8));
vsum = _mm_add_epi32(vsum, _mm_srli_si128(vsum, 4));
result = _mm_cvtsi128_si32(vsum);
// handle any residual bytes ( < 16)
if (i < length)
{
result += fast_compare_ref(&s[i], &t[i], length - i);
}
return result;
}
// test harness
int main(void)
{
const int n = 1000000;
char *s = malloc(n);
char *t = malloc(n);
int i, result_ref, result;
srand(time(NULL));
for (i = 0; i < n; ++i)
{
s[i] = rand();
t[i] = rand();
}
result_ref = fast_compare_ref(s, t, n);
result = fast_compare(s, t, n);
printf("result_ref = %d, result = %d\n", result_ref, result);;
return 0;
}
Compile and run the above test harness:
$ gcc -Wall -O3 -msse3 fast_compare.c -o fast_compare
$ ./fast_compare
result_ref = 3955, result = 3955
$ ./fast_compare
result_ref = 3947, result = 3947
$ ./fast_compare
result_ref = 3945, result = 3945
Note that there is one possibly non-obvious trick in the above SSE code where we use _mm_madd_epi16
to unpack and accumulate 16 bit 0
/-1
values to 32 bit partial sums. We take advantage of the fact that -1*-1 = 1
(and 0*0 = 0
of course) - we're not really doing a multiply here, just unpacking and summing in one instruction.
UPDATE: as noted in the comments below, this solution is not optimal - I just took a fairly optimal 16 bit solution and added 8 bit to 16 bit unpacking to make it work for 8 bit data. However for 8 bit data there are more efficient methods, e.g. using psadbw
/_mm_sad_epu8
. I'll leave this answer here for posterity, and for anyone who might want to do this kind of thing with 16 bit data, but really one of the other answers which doesn't require unpacking the input data should be the accepted answer.