suffix-array.md

Суффиксный массив

Суффиксный массив, автомат и дерево обобщённо называют суффиксными структурами данных. Они применяются в множестве различных задач, встречающихся как на олимпиадах, так и на практике.

Суффиксные структуры часто (но не всегда) взаимозаменяемые, и более того, конвертируются друг в друга за линейное время. Суффиксный массив — самый простой из них, поэтому мы с него и начнём.

«*Паблик с тупыми шутками про проганье*»

Мотивация

Суффиксным массивом строки $s$ называется перестановка индексов начал её суффиксов, которая задаёт их порядок в порядке лексикографической сортировки. Иными словами, чтобы его построить, нужно выполнить сортировку всех суффиксов заданной строки.

Как это использовать. Пусть вы решили основать ООО «Ещё Один Поисковик», и чтобы получить финансирование, вы хотите сделат хоть что-то минимально работающие — просто научиться искать по ключевому слову документы, включающие его, а также позиции их вхождения (в 90-е это был бы уже довольно сильный MVP). Простыми алгоритмами (полиномиальными хэшами, z- и префикс-функцией и даже Ахо-Корасиком) это сделать быстро нельзя, а суффиксными структурами — можно.

В случае с суффиксным массивом можно сделать следующее: найти бинарным поиском первый суффикс в суффиксном массиве, который меньше искомого слова, а также последний, который меньше. Все суффиксы между этими двумя будут включать искомую строку как префикс.

Работать такой алгоритм будет за $O(|t| \log |s|)$, и позже это можно будет оптимизировать до $O(|t| + \log |s|)$, что является одним из самых оптимальных алгоритмов поиска.

Теперь научимся его строить.

Построение за $O(n \log n)$

Для удобства мы допишем в конец строки какой-нибудь символ, который лексикографически меньше любого другого, и будем выполнять сортировку не суффиксов, а циклических сдвигов. Строки мы, соответственно, тоже будем рассматривать циклические. Для стандартной ASCII обычно выбирают либо «$» либо «#», либо просто нулевой символ, если это c-string (в языке Си все строки и так заканчиваются нулевым символом). Легко убедиться, что сортировка таких циклических сдвигов эквивалентна сортировке суффиксов — можно просто убрать всё, что идёт после доллара.

Мы могли бы просто взять перестановку от $0$ до $n$, написать компаратор, который сравнивает соответствующие суффиксы, и скормить это в std::sort, что будет работать за $O(n^2 \log n)$, потому что внутреннее сравнение работает за $O(n)$. Однако, если сравнивать суффиксы хэшами, то уже тут можно получить $O(n \log^2 n)$. Но это не самый быстрый и удобный алгоритм.

Наш алгоритм будет состоять из $\lceil \log n \rceil$ этапов. На $k$-том этапе мы будем рассматривать циклические подстроки длины $2^k$. На последнем этапе мы отсортируем строки длины $\geq n$ (это легально — они ведь циклические), и мы получим нужный суффиксный массив.

Заметим, что, в отличие сортировки суффиксов, сортировка подстрок не всегда однозначная — они могут быть одинаковыми. Поэтому на каждой фазе алгоритм помимо перестановки $p$ индексов циклических подстрок мы будем поддерживать для каждой циклической подстроки, начинающейся в позиции $i$ с длиной 2^k, номер $c_i$ класса эквивалентности, которому эта подстрока принадлежит (давать их будем таким образом, чтобы они сохраняли порядок: меньшим подстрокам соответствуют меньшие $c_i$). Количество классов эквивалентности будем хранить в переменной cls (изначально она равна количеству различных символов).

Пример: $s = aaba$. Этапов будет 3: для подстрок длины 1, 2 и 4.

$$ p_0 = (0, 1, 3, 2) ;;; c_0 = (0, 0, 1, 0) \ p_1 = (0, 3, 1, 2) ;;; c_1 = (0, 1, 2, 0) \ p_2 = (3, 0, 1, 2) ;;; c_2 = (1, 2, 3, 0) $$

Как настоящие программисты, мы нумеруем этапы с нуля, поэтому на нулевом этапе мы отсортируем строки длины $2^0 = 1$, то есть просто символы. Это мы сделаем сортировкой подсчётом.

Следующие этапы нужно проводить, используя информацию с предыдущих. Можем снова создать перестановку и применить к ней std::sort со своим компаратором.

Как быстро сравнить две подстроки? Мы можем использовать $c_i$ — каждой строке длины $2^k$ сопоставить биграмму (строку из двух символов), а именно строка $s[i..i+2^k-1]$ с точки зрения сортировки будет эквивалентна паре $(c_i, c_{i+2^{k-1}})$. Соответственно, можно просто написать компаратор, который смотрит только на эти пары, и работает за $O(1)$. Однако, это всё ещё будет работать за $O(n \log^2 n)$, потому что каждый этап будет работать за $O(n \log n$).

Примечание. Зачастую этот способ построения работает быстро, поскольку имеет небольшую константу.

Оптимизация до $O(n \log n)$. Сортировка слиянием оптимальна только тогда, когда мы ничего не знаем про сортируемые значения — в нашем случае это не так. Воспользуемся цифровой сортировкой — сначала отсортируем биграммы по второму символу, а потом по первому. Все вторые элементы уже упорядочены (массив $p$ с предыдущего этапа). Теперь, чтобы упорядочить их по первому, нам надо просто от каждого элемента в $p$ отнять $2^{k-1}$. Таким образом, можно проводить этап за $O(n)$.

// строка -- это последовательность чисел от 1 до размера алфавита
vector<int> suffix_array (vector<int> &s) {
    s.push_back(0);  // добавляем нулевой символ в конец строки
    int n = (int) s.size(),
        cnt = 0,  // вспомогательная переменная: счётчик для сортировки 
        cls = 0;  // количество классов эквивалентности
    vector<int> c(n), p(n);
    
    map< int, vector<int> > t;
    for (int i = 0; i < n; i++)
        t[s[i]].push_back(i);
    
    // «нулевой» этап
    for (auto &x : t) {
        for (int u : x.second)
            c[u] = cls, p[cnt++] = u;
        cls++;
    }
    
    // пока все суффиксы не стали уникальными
    for (int l = 1; cls < n; l++) {
        vector< vector<int> > a(cls);  // массив для сортировки подсчётом
        vector<int> _c(n);  // новые классы эквивалентности
        int d = (1<<l)/2;
        int _cls = cnt = 0;  // новое количество классов
        
        for (int i = 0; i < n; i++) {
            int k = (p[i]-d+n)%n;
            a[c[k]].push_back(k);
        }
        
        for (int i = 0; i < cls; i++) {
            for (size_t j = 0; j < a[i].size(); j++) {
                // если суффикс начинает новый класс эквивалентности
                if (j == 0 || c[(a[i][j]+d)%n] != c[(a[i][j-1]+d)%n])
                    _cls++;
                _c[a[i][j]] = _cls-1;
                p[cnt++] = a[i][j];
            }
        }
        
        c = _c;
        cls = _cls;
    }
    
    return vector<int>(p.begin()+1, p.end());
}

TODO: переписать это

Наибольшие общие префиксы

Для многих применений очень часто требуется искать длину наибольшего общего префикса (англ — largest common prefix) для различных суффиксов строки. Например, просто чтобы сравнить две подстроки, мы можем найти общий префикс соответствующих суффиксов и посмотреть на следующий символ — так же, как мы делали с хэшами.

Пусть мы знаем $lcp$ для всех соседей в суффиксном массиве. Тогда для того, чтобы найти lcp(s[i:], s[j:]), надо найти $\min_{k=c[i]}^{c[j-1]}$lcp(s[p[k]:], s[p[k + 1]:]) (для доказательства можно представить массив $lcp$ в виде гистограммы). Тогда есть мотивация посчитать массив $lcp$, в котором окажутся наибольшие общие префиксы соседних суффиксов, а после считать минимумы на отрезках в этом массиве (к примеру, с помощью разреженных таблиц).

Осталось придумать способ быстро посчитать массив $lcp$. Можно воспользоваться идеей из построения суффиксного массива за $O(n \log^2 n)$ — с помощью хэшей и бинпоиска находить $lcp$ для каждой пары соседей. Такой метод работает за $O(n \log n)$, но является не самым удобным и популярным.

Алгоритм Касаи, Аримуры, Арикавы, Ли, Парка. В простонародье называется как угодно, но не исходным способом (алгоритм Касаи, алгоритм пяти корейцев, итд.). Используется для подсчета $lcp$ за $O(n)$. Автору алгоритм кажется чем-то похожим на z-функцию по своей идее. Пусть мы уже посчитали $lcp[i]$ (считаем, что $i$ — позиция суффикса в уже построенном суфмассе). Заметим, что $lcp[c[p[i] + 1]] \ge lcp[i] - 1$.

Доказательство — если для $p[i]$-го суффикса была общая часть длины $lcp[i]$ с $p[i + 1]$-ым суффиксом, то у $p[i] + 1$-го суффикса гарантированно будет общая часть длины $lcp[i] - 1$ с $p[i + 1] + 1$-ым суффиксом. Заметим, что этот суффикс уже не обязательно будет соседом $p[i] + 1$-го в суффиксном массиве, но из свойства $lcp$ двух произвольных суффиксов получим, что $lcp[c[p[i] + 1]] \ge lcp[i] - 1$. Вычислять точное значение $lcp$ мы будем наивным образом — увеличиваем, пока увеличивается, а вот брать стартовое значение будем, опираясь на это неравенство. В таком случае правильный порядок вычисления массива $lcp$ — от длинных суффиксов к коротким.

Для того, чтобы оценить сложность алгоритма, посчитаем, сколько раз мы сделаем наивное прибавление. $lcp[i] \le n$, а также мы сделаем не более $n$ вычитаний единицы при переходах. Тогда суммарное число прибавлений — $O(n)$. Следовательно, и сам алгоритм будет работать за $O(n)$.

Код алгоритма:

    vector<int> calc_lcp(vector<int> &val, vector<int> &c, vector<int> &p) {
        int n = val.size();
        int current_lcp = 0;
        vector<int> lcp(n);
        for (int i = 0; i < n; i++) {
            if (c[i] == n - 1)
                continue;
            int nxt = p[c[i] + 1];
            while (max(i, nxt) + current_lcp < n && val[i + current_lcp] == val[nxt + current_lcp])
                current_lcp++;
            lcp[c[i]] = current_lcp;
            current_lcp = max(0, current_lcp - 1);
        }
        return lcp;
    }