code

Crochemore-Perrin algorithm using SSE instructions

The Crochemore-Perrin algorithm using SSE instructions algorithm
ssecp.c
#include "include/log2.h"
#include <assert.h>
#include <nmmintrin.h>


const int max_needle = 8;

int search_rawsse(unsigned char *x, int m, unsigned char *y, int n) 
{
    if (m > max_needle)
return -1;

    int occurences = 0;
    __m128i needle_reg = _mm_loadu_si128((__m128i *)x);

    int step = 16 - m + 1;     unsigned int result_mask = (1 << step) - 1;

    int number_of_steps = (n + step - 16) / step;
    int steps_size = number_of_steps * step;
    unsigned char *steps_end = y + steps_size;

    while (y != steps_end) { 
        __m128i haystack_reg = _mm_loadu_si128((__m128i *)y);
        __m128i mask_reg = _mm_cmpestrm(needle_reg, m, haystack_reg, 16, SIDD_UBYTE_OPS | SIDD_CMP_EQUAL_ORDERED);
        occurences += _mm_popcnt_u32(*(int *)&mask_reg & result_mask);
        y += step;
    }

    n -= steps_size; 
    if (n >= m) {
        __m128i haystack_reg = _mm_loadu_si128((__m128i *)y);
        __m128i mask_reg = _mm_cmpestrm(needle_reg, m, haystack_reg, n, SIDD_UBYTE_OPS | SIDD_CMP_EQUAL_ORDERED);
        step = n - m + 1;
        result_mask = (1 << step) - 1;
        occurences += _mm_popcnt_u32(*(int *)&mask_reg & result_mask);
    }

    return occurences;
}

int maxSuf(char *x, int m, int *p) {
   int ms, j, k;
   char a, b;

   ms = -1;
   j = 0;
   k = *p = 1;
   while (j + k < m) {
      a = x[j + k];
      b = x[ms + k];
      if (a < b) {
         j += k;
         k = 1;
         *p = j - ms;
      }
      else
         if (a == b)
            if (k != *p)
               ++k;
            else {
               j += *p;
               k = 1;
            }
         else { /* a > b */
            ms = j;
            j = ms + 1;
            k = *p = 1;
         }
   }
   return(ms);
}
 
int maxSufTilde(char *x, int m, int *p) {
   int ms, j, k;
   char a, b;

   ms = -1;
   j = 0;
   k = *p = 1;
   while (j + k < m) {
      a = x[j + k];
      b = x[ms + k];
      if (a > b) {
         j += k;
         k = 1;
         *p = j - ms;
      }
      else
         if (a == b)
            if (k != *p)
               ++k;
            else {
               j += *p;
               k = 1;
            }
         else { /* a < b */
            ms = j;
            j = ms + 1;
            k = *p = 1;
         }
   }
   return(ms);
}

void preKmp(unsigned char *x, int m, int kmpNext[]) {
   int i, j;
   i = 0;
   j = kmpNext[0] = -1;
   while (i < m) {
      while (j > -1 && x[i] != x[j])
         j = kmpNext[j];
      i++;
      j++;
      if (i<m && x[i] == x[j])
         kmpNext[i] = kmpNext[j];
      else
         kmpNext[i] = j;
   }
}

void  compute(unsigned char* x, int m, int *mu, int *pi)
{
   int i, j, ell, p, per, q, kmpNext[XSIZE];

   /* Preprocessing */
   i = maxSuf(x, m, &p);
   j = maxSufTilde(x, m, &q);
   if (i > j) {
      ell = i;
      per = p;
   } else {
      ell = j;
      per = q;
   }

   *mu = ell+1;

   if (!memcmp(x,x + per, m-per)) {
       *pi = per;
   } else {
       int kmpNext[XSIZE];
       preKmp(x, m, kmpNext);
       *pi = m-kmpNext[m];
   }
}


int search(unsigned char *x, int m, unsigned char *y, int n) 
{
    if (m <= max_needle) 
return search_rawsse(x, m, y, n);

    int mu, pi, count = 0;
    compute(x, m, &mu, &pi);

    assert( (mu > 0 && mu < pi && mu < m  && pi > 0 && pi <= m) );

    __m128i needle_reg = _mm_loadu_si128((__m128i *)&x[mu]);
    int needle_length = m - mu;
    if (needle_length > 16)
needle_length = 16;

loop:  
    while (n >= m) {
        __m128i haystack_reg = _mm_loadu_si128((__m128i *)&y[mu]);

int haystack_length = n - mu;
if (haystack_length > 16)
haystack_length = 16; 

        int idx = _mm_cmpestri(needle_reg, needle_length, haystack_reg, haystack_length, SIDD_UBYTE_OPS | SIDD_CMP_EQUAL_ORDERED);
if (!_mm_cmpestrc(needle_reg, needle_length, haystack_reg, haystack_length, SIDD_UBYTE_OPS | SIDD_CMP_EQUAL_ORDERED)) {
y += haystack_length;
n -= haystack_length;

goto loop; 
}

y += idx;
n -= idx;
if (n < m) 
break; 

int head0 = 0, head;

if (needle_length <= haystack_length - idx) {
head = mu + needle_length; 
} else {
head = mu + haystack_length - idx; 
}

loop1:
while (head < m) {
__m128i b0 = _mm_loadu_si128((__m128i *)&x[head]), 
b1 = _mm_loadu_si128((__m128i *)&y[head]);

unsigned int b = m - head;

unsigned int idx = _mm_cmpestri(b0, b, b1, b, SIDD_UBYTE_OPS | SIDD_CMP_EQUAL_EACH | SIDD_NEGATIVE_POLARITY);
if (_mm_cmpestrc(b0, b, b1, b, SIDD_UBYTE_OPS | SIDD_CMP_EQUAL_EACH | SIDD_NEGATIVE_POLARITY)) {
             
head += idx + 1;
int step = head - mu;
y += step;
n -= step;

goto loop; 
}

head += idx;
}

loop2:
while (head0 < mu) {
__m128i b0 = _mm_loadu_si128((__m128i *)&x[head0]), 
b1 = _mm_loadu_si128((__m128i *)&y[head0]);

unsigned int b = mu - head0;

unsigned int idx = _mm_cmpestri(b0, b, b1, b, SIDD_UBYTE_OPS | SIDD_CMP_EQUAL_EACH | SIDD_NEGATIVE_POLARITY);
if (_mm_cmpestrc(b0, b, b1, b, SIDD_UBYTE_OPS | SIDD_CMP_EQUAL_EACH | SIDD_NEGATIVE_POLARITY)) {
goto loop3;
}

head0 += idx;
}

count++;

loop3:
y += pi;
n -= pi;
head = m - pi;
if (head < mu) { 
head0 = head;
head = mu;
}

goto loop1;
    }

    return ++count;
}