Оптимизация скорости исполнения кода

Данная статья о том, как классические алгоритмы позволяют сделать наш код быстрее. Мы рассмотрим алгоритм поиска нулевого бита и то как он нам может помочь повысить эффективность алгоритма поиска символа из диапазона (например найти первое вхождение символа из диапазона 0-9 в строке).
Т.е. просто сферический конь в вакууме довольно распространенная ситуация при обработке текста, когда необходимо найти положение первой цифры в произвольном тексте. То, с чего многие начинают учить регулярные выражения.
В статье описывается использование gcc. При желании все можно повторить под другие компиляторы и языки с небольшими переделками. Осторожно, под катом есть ассемблер!

Для решения этой задачи мы можем воспользоваться такими путями:
1) Воспользоваться библиотекой регулярных выражений.
2) Написать функцию для поиска самостоятельно.
В случае если строки маленькие, вопрос об оптимизации не возникает. Данная задача становиться интересной когда длинна строк увеличивается, либо имеет место большое количество проверок.
Время исполнения будем замерять при помощи функции clock() из подключаемого файла <time.h>.

Используем REGEX

Рассмотрим способ с использованием регулярных выражений. Один из первых методов который приходит в голову — воспользоваться библиотекой REGEX. Регулярное выражение будет выглядеть так:
"[[:digit:]]"
Код для использования этой библиотеки может быть таким:

clock_t tbefore;
    long i;    
    double eplased;    
    int pmatch[1];
    regex_t regex;
    int reti;
….
/*REGEXP usage*/
   tbefore = clock();
   reti = regcomp( & regex, "[[:digit:]]", 0);   
   for (i=1;i<100000;i++){
        reti = regexec(& regex, test_string, 1, pmatch, 0);	
   }   
   eplased=clock()-tbefore;
   printf("Function strchr used %.3f seconds %d\n", eplased/CLOCKS_PER_SEC, pmatch[0]);    
   regfree( & regex);
/**************/

Скорость исполнения: 1.74 сек.
В данном коде я намеренно убрал проверки на успешность компиляции и поиска, т. е. использовать данный код для решения реальных задач нельзя. Но для проверки скорости в «тепличных» условиях — самое то. Кроме того, я замеряю время исполнения с учетом компиляции регулярного выражения и освобождения памяти после его использования.

Своя функция для поиска

Попробуем воспользоваться традиционным подходом — в лоб. Функция поиска
будет выглядеть таким образом (Алгоритм 1):

int search1(char *str){
  int l=0;
  int i;
  
  l = strlen(str);
  for (i=0;i<l;i++){
    if ((str[i] >= '0') && (str[i] <= '9')){
      return (i);
    }
  }  
  
  return (-1);
}

Скорость исполнения: 7.19 сек.
И результаты замеров показывают более чем в четыре раза меньшую скорость исполнения! Чего и следовало ожидать.

Оптимизированный алгоритм

Попробуем оптимизировать. Когда-то в руки мне попалась книга «Алгоритмические трюки для программистов» Г. Уоррена. И в этой книге было описание алгоритмов поиска нулевого символа в строке. Эти алгоритмы используются например в языке С для поиска конца строки. Также в одном из абзацев предложено обобщение одного из таких алгоритмов для описка диапазона.
Идея алгоритма заключается в том, что для поиска индекса символа в строке перебирают не по по дному символу, а сразу по четыре, считывая из массива двойные слова(dword = 4 byte). Т.е. В терминах языка С мы из массива символов (char) считываем длинные целые (unsigned).
После этого этого со словом проводятся арифметические действия. Для диапазона от нуля до девяти это выглядит таким образом:

int zbytel3(unsigned z) {
   unsigned y;
   x = z ^ 0x30303030;  
                         				// Original byte: 00 80 other
   y = (x & 0x7F7F7F7F) + 0x76767676; // 7F 7F 1xxxxxxx
   y = ~(y | x | 0x7F7F7F7F);         // 80 00 00000000
                                      // These steps map:
   if (y == 0) return 4;              // 00000000 ==> 4,
   else if (y > 0x0000FFFF)           // 80xxxxxx ==> 0,
      return (y >> 31) ^ 1;           // 0080xxxx ==> 1,
   else                               // 000080xx ==> 2,
      return (y >> 15) ^ 3;           // 00000080 ==> 3.
}
    

т. е. x = z ^ 0x30303030 — операция исключающее или, которая превращает минимальное число диапазона „0“ — 0x30 в 0. 0x76 = 0x7F-9, где 9=n-1, где n — количество символов в диапазоне.
Далее простая байтовая арифметика, в результате которой можно получить число, которое обрабатыватеся в ветвлении, и в результате мы получаем цифры от 1 до 4. Т.е. если символ не найден — то 4.
Теперь наша функция примет такой вид:

int search2( char *str){
  int n, j = 0;  
  unsigned x, y, *r;
  const char *c_ptr;
  //Проверяем посимвольно первые байты, пока не выровняем.
   for (c_ptr=str; ((unsigned long int)c_ptr & (sizeof(x) - 1)) != 0; ++c_ptr)
    if ((*c_ptr >= '0') && (*c_ptr <= '9'))
      return c_ptr - str;
  
  r = (unsigned *) c_ptr;
  j = c_ptr - str ;  
  while (1){            
      x = *r ^ 0x30303030;                                           	    
      y = (x & 0x7F7F7F7F) + 0x76767676;
      y = ~(y | x | 0x7F7F7F7F);         
      
      // These steps map:
      if (y == 0) n = 4;              // 00000000 ==> 4,
      else if (y > 0x0000FFFF)        // 80xxxxxx ==> 0,
      n= (y >> 31) ^ 1;               // 0080xxxx ==> 1,
      else                            // 000080xx ==> 2,
      n= (y >> 15) ^ 3;               // 00000080 ==> 3.            
      j=j+n ;
      if (n<4) {j =j +3 -2*n; break;}      
      r++;
  }        
  return (j);
}

Скорость исполнения: 1.71 сек.
Результаты замеров показывают скорость близкую к regexp. Можем ли мы сделать быстрее? Да!
Воспользуемся ассемблерными вставками. Сила описанного метода в том, что он позволяет уменьшить количество тактов которые тратит процессор на выполнение за счет применения элементарных арифметических операций. Программирование на языке высокого уровня не до конца позволяет реализовать возможности этого алгоритма.
А потому, теперь наша функция примет такой вид:

int search3( char *str){
  int n, j = 0;
  unsigned x, y, *r;
  const char *c_ptr;    
  for (c_ptr=str; ((unsigned long int)c_ptr & (sizeof(x) - 1)) != 0; ++c_ptr)
    if ((*c_ptr >= '0') && (*c_ptr <= '9'))
      return c_ptr - str;
   
  r = (unsigned *) c_ptr;
  j = c_ptr - str ;
  while (1){
      __asm__(
       "movl $5, %%edx\n"
       "movl (%%ebx), %%eax\n"
       "xor $0x30303030, %%eax\n" //
       "and $0x7F7F7F7F, %%eax\n"
       "add $0x76767676, %%eax\n"
       "movl %%eax, %%ecx\n"
       "or %%eax, %%ecx\n"
       "or $0x7F7F7F7F, %%ecx\n"
       "notl %%ecx\n"                   
      :"=c"(y)
      :"b"( r )
      );
      
       // These steps map:
      if (y == 0) n = 4;              // 00000000 ==> 4,
      else if (y > 0x0000FFFF)        // 80xxxxxx ==> 0,
      n= (y >> 31) ^ 1;               // 0080xxxx ==> 1,
      else                            // 000080xx ==> 2,
      n= (y >> 15) ^ 3;               // 00000080 ==> 3.            
      j=j+n ;
      if (n<4) {j =j +3 -2*n; break;}      
      r++;
  }        
  return (j);
}

Скорость исполнения: 1.23 сек.
Я заменил ассемблерной вставкой только арифметические операции. Ветвление тоже можно заменить ассемблерным кодом, но я не уверен в значительном выигрыше, но несколько ветвлений с метками значительно усложнят процесс отладки.
Выигрыш присутствует, но не в разы по сравнению с regexp. В принципе, оптимизируя ассемблерными вставками можно еще улучшать код, но это все влечет усложнение процесса отладки. Если оптимизировать ветвления и цикл, то можно добиться увеличения скорости еще раза в полтора. Но отлаживать код будет значительно труднее. Еще можно оптимизировать этот алгоритм для 64 бит вместо 32, попытаться распараллелить процесс.

Итого, получаем:

Regex	1.74
Алгоритм 1	7.19
Алгоритм 2	1.71
Алгоритм 2+asm	1.23

Вывод напрашивается сам: Изобретение велосипедов не приветствуется, и не всегда выигрыш от ассемблерной оптимизации может быть столь значим.

Ссылка на архив с исходниками:
webfile.ru/5706721

UPD:
Оптимизация под архитектуру процессора была указана, добавляю параметр -O 2.
Также добавляю суммирование внутрь тестовых функций, чтобы как-либо использовать результат.
Результат:

Regex	4.56
Алгоритм 1	2.99
Алгоритм 2	1
Алгоритм 2+asm	1.7*

*С ключем -O2 правильно работать перестал, __volatile_ не помогло, результат с ключем -O1 но с добавленным суммированием.

UPD 2:

Пользователем Maratyszcza (спасибо ему!) предложен еще более быстрый вариант с использованием SIMD-расширений.

#include <intrin.h>

const char* find_digit(const char* str) {
    static const __m128i str_mask[16] = {
        _mm_set_epi32(0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF),
        _mm_set_epi32(0x00FFFFFF, 0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF),
        _mm_set_epi32(0x0000FFFF, 0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF),
        _mm_set_epi32(0x000000FF, 0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF),
        _mm_set_epi32(0x00000000, 0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF),
        _mm_set_epi32(0x00000000, 0x00FFFFFF, 0xFFFFFFFF, 0xFFFFFFFF),
        _mm_set_epi32(0x00000000, 0x0000FFFF, 0xFFFFFFFF, 0xFFFFFFFF),
        _mm_set_epi32(0x00000000, 0x000000FF, 0xFFFFFFFF, 0xFFFFFFFF),
        _mm_set_epi32(0x00000000, 0x00000000, 0xFFFFFFFF, 0xFFFFFFFF),
        _mm_set_epi32(0x00000000, 0x00000000, 0x00FFFFFF, 0xFFFFFFFF),
        _mm_set_epi32(0x00000000, 0x00000000, 0x0000FFFF, 0xFFFFFFFF),
        _mm_set_epi32(0x00000000, 0x00000000, 0x000000FF, 0xFFFFFFFF),
        _mm_set_epi32(0x00000000, 0x00000000, 0x00000000, 0xFFFFFFFF),
        _mm_set_epi32(0x00000000, 0x00000000, 0x00000000, 0x00FFFFFF),
        _mm_set_epi32(0x00000000, 0x00000000, 0x00000000, 0x0000FFFF),
        _mm_set_epi32(0x00000000, 0x00000000, 0x00000000, 0x000000FF)
    };
    static const __m128i str_offset = _mm_set1_epi8(127 - '9');
    static const __m128i str_threshold = _mm_set1_epi8(127 - 10);
    const size_t str_misalignment = ((size_t)str) & ((size_t)15);
    const __m128i* str_aligned = (const __m128i*)(((size_t)str) - str_misalignment);
    __m128i str_vector = _mm_load_si128(str_aligned);
    str_vector = _mm_and_si128(str_vector, str_mask[str_misalignment]);
    str_vector = _mm_add_epi8(str_vector, str_offset);
    int str_bitmask = _mm_movemask_epi8(_mm_cmpgt_epi8(str_vector, str_threshold));
    unsigned long index;
    _BitScanForward(&index, str_bitmask);
    while (str_bitmask == 0) {
        str_aligned += 1;
        str_vector = _mm_load_si128(str_aligned);
        str_vector = _mm_add_epi8(str_vector, str_offset);
        str_bitmask = _mm_movemask_epi8(_mm_cmpgt_epi8(str_vector, str_threshold));
        _BitScanForward(&index, str_bitmask);
    }
    return ((const char*)str_aligned) + index;
}

Пока это самый быстрый вариант.