Двоичный алгоритм поиска подстроки

Двоичный алгоритм поиска подстроки (также bitap algorithm, shift-or algorithm) — алгоритм поиска подстроки, использующий тот факт, что в современных компьютерах битовый сдвиг и побитовое ИЛИ являются атомарными операциями. По сути, это примитивный алгоритм поиска с небольшой оптимизацией, благодаря которой за одну операцию производится до 32 сравнений одновременно (или до 64, в зависимости от разрядности машины). Легко переделывается на приблизительный поиск.

Вычислительная сложность — O(|needle|·|haystack|) операций с крайне малой константой. На предварительную обработку идёт O(|needle|·|Σ|) операций и памяти, где Σ — алфавит. Если же длина шаблона поиска (в символах) не превышает длину машинного слова (в битах), сложность получается O(|haystack|) и O(|needle|+|Σ|) соответственно.

Предполагается, что алгоритм разработал в 1964 году венгр Балинт Дёмёльки. Но популярность он приобрёл в 1990-е годы благодаря работам чилийца Рикардо Баэса-Ятеса и уругвайца Гастона Гоннета. На двоичном алгоритме основана известная утилита приблизительного поиска agrep.

Идея

Как всегда, считаем, что needle («иголка») — строка, которую мы ищем, а haystack («стог сена») — строка, в которой ищем. Длины этих строк обозначим через n и H. Символы в строке пронумерованы с 1, как в Паскале.

Построим такую матрицу n×H: если i-префикс needle совпадает с haystack в позициях j−i+1…j, ставим в матрице 1 в позиции (i, j). Иначе — ноль.

      haystack
    | Г Ц А Т  Ц Г Ц А  Г А Г А  Г Т А Т  А Ц А Г  Т А Ц Г
  --------------------------------------------------------
n Г | 1 0 0 0  0 1 0 0  1 0 1 0  1 0 0 0  0 0 0 1  0 0 0 1
e Ц | 0 1 0 0  0 0 1 0  0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0
e А | 0 0 1 0  0 0 0 1  0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0
d Г | 0 0 0 0  0 0 0 0  1 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0
l А | 0 0 0 0  0 0 0 0  0 1 0 0  0 0 0 0  0 0 0 0  0 0 0 0
e Г | 0 0 0 0  0 0 0 0  0 0 1 0  0 0 0 0  0 0 0 0  0 0 0 0
  А | 0 0 0 0  0 0 0 0  0 0 0 1  0 0 0 0  0 0 0 0  0 0 0 0
  Г | 0 0 0 0  0 0 0 0  0 0 0 0  1 0 0 0  0 0 0 0  0 0 0 0

Шаблон найден, если найдена единица в последней строке этой матрицы. Матрицу будем вычислять динамически, по столбцам. Для этого зададимся вопросом: как превратить (j−1)-й столбец матрицы в j-й?

R[1,j] = 1, если needle[1] = haystack[j].

R[2,j] = 1, если needle[2] = haystack[j] и R[1,j−1] = 1.

R[3,j] = 1, если needle[3] = haystack[j] и R[2,j−1] = 1.

R[4,j] = 1, если needle[4] = haystack[j] и R[3,j−1] = 1.

…

R[n,j] = 1, если needle[n] = haystack[j] и R[n−1,j−1] = 1.

Запишем первую формулу как

R[1,j] = 1, если needle[1] = haystack[j] и ИСТИНА.

Если объединить всё это в один двоичный кортеж R{j} длины n, получится такая формула:

R{j} = ((R{j−1} << 1) | 00…001) & comparison_vector[haystack[j]].

Операция << — двоичный сдвиг влево, операция & — побитовое И, | — побитовое ИЛИ. Вектор comparison_vector строится так: в элементе $a\in \Sigma$ на i-й позиции стоит 1, если в needle на этой позиции стоит символ a.

      алфавит
    | А Г Т Ц
  -----------
n Г | 0 1 0 0
e Ц | 0 0 0 1
e А | 1 0 0 0
d Г | 0 1 0 0
l А | 1 0 0 0
e Г | 0 1 0 0
  А | 1 0 0 0
  Г | 0 1 0 0

В этом и заключается предварительная обработка подстроки needle.

Поскольку на реальных ЭВМ при двоичном сдвиге на освободившееся место приходит 0, часто единицу и ноль меняют ролями. Формула получается

T{j} = (T{j−1} << 1) | comparison_vector_neg[haystack[j]].

Или, на языке программирования:

T = (T << 1) | comparison_vector_neg[haystack[j]];

Отсюда название «Shift-Or».

Исходный текст

Си

Примечание. Работаем с «перевёрнутыми» битами (0 — совпадает, 1 — нет).

// Размер алфавита
#define ASIZE  256
// Размер машинного слова
#define WORD (sizeof(unsigned int)*8)

int preSo(char *x, int m, unsigned int S[]) { 
  unsigned int j, lim; 
  int i; 
  for (i = 0; i < ASIZE; ++i) 
    S[i] = ~0u; 
  for (lim = i = 0, j = 1; i < m; ++i, j <<= 1) { 
    S[x[i]] &= ~j; 
    lim |= j; 
  } 
  lim = ~(lim>>1); 
  return(lim); 
}

/*
  x - needle, длина m
  y - haystack, длина n
*/
void SO(char *x, int m, char *y, int n) { 
  unsigned int lim, state; 
  unsigned int S[ASIZE]; 
  int j; 
  if (m > WORD) 
    error("SO: шаблон поиска слишком длинный");

  /* Предварительная обработка */ 
  lim = preSo(x, m, S);

  /* Поиск */
  for (state = ~0, j = 0; j < n; ++j) { 
    state = (state<<1) | S[y[j]]; 
    if (state < lim) 
      OUTPUT(j - m + 1); 
  } 
}

Версия для приблизительного поиска

Для простоты работаем с «перевёрнутыми» битами (1 — не совпадает, 0 — совпадает). Предполагается, что needle содержит не более m ошибок в метрике Левенштейна. Нам нужна m+1 копия кортежа T.

q{j, k} = (T{j−1, k} << 1) | comparison_vector_neg[haystack[j]];

Ins{j, k} = q{j, k} & T{j−1, k−1};

Del{j, k} = q{j, k} & (T{j, k−1} << 1);

Repl{j, k} = q{j, k} & (T{j−1, k−1} << 1);

T{j, k} = Ins{j, k} & Del{j, k} & Repl{j, k},

где j=1…H, k=0…m. Вектор T(*, −1) состоит из одних единиц. Поиск удачен, если хотя бы в одном из векторов T в строке n окажется ноль. Соответствующее k — количество ошибок.

Сравнение с другими алгоритмами

Преимущества

Очень быстр на подстроках, длина которой (в символах) не превышает длину машинного слова (в битах). Нет никакой регрессии на «плохих» данных.

Алгоритм легко переделывается на приблизительный поиск в метрике Хэмминга и даже в метрике Левенштейна, а также на поиск одной из нескольких строк.

Недостатки

При точном поиске не имеет особых преимуществ над другими алгоритмами (например, Бойера—Мура).

Непонятно, что делать на громадных алфавитах (например, Юникод). Например, можно разбить один «длинный» символ на несколько, при этом в позициях, некратных длине символа, в R входит не единица, а ноль. Или можно воспользоваться какой-то структурой данных, которая хранила бы «разреженный» массив.

При приблизительном поиске алгоритм указывает конец места, где совпало. Начало этого места найти сложнее, требуется O(m·n) дополнительной памяти.

Версия для длинных подстрок программируется несколько сложнее. Понятия «флаг переноса» в языках высокого уровня нет, поэтому приходится либо рассчитывать на удачную оптимизацию, либо оптимизировать код вручную на ассемблере.