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

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

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

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

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

Мотивация

Суффиксным массивом строки \(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;
    }