Сумма по подмножествам и xor-and-or-свертки
TL;DR: SOS-DP — это алгоритм, который за $O(2^k \cdot k)$ находит для каждой из $2^k$ масок сумму элементов изначального массива по всем подмаскам. Считается при помощи динамики, перебирая биты и постепенно пытаясь их менять. or-xor-and-свертки — это обобщения свертки двух последовательностей через преобразование Фурье, но вместо $c_i = \sum_{j, k : j + k = i} a_j \cdot b_k$ мы сворачиваем $c_i = \sum_{j, k : j \diamondsuit k = i} a_j \cdot b_k$, где $\diamondsuit$ — это какая-то битовая операция (or, xor, and). Все свертки получаются применением какого-то преобразования к исходным последовательностям, поэлементным их перемножением и применением обратного преобразования. Для or-свертки это преобразование — сумма по подмножествам, для and-свертки — сумма по надмножествам, для xor-свертки — преобразование Уолша — Адамара. Во всех случаях, чтобы получить обратное преобразование, необходимо просто обратить все действия, которые делает прямое преобразование, в обратном порядке.
Сумма по подмножествам
Определение: Пусть дан массив $a$. Назовем массивом сумм по подмножествам (Sum Over Subsets или SOS-DP) массива $a$ такой массив $b$, что
$$ b_i = \sum_{j : j \subset i} a_j$$
При этом условие $j \subset i$ в данном случае воспринимается как содержание подмножеств единичных битов. То есть если бы мы воспринимали числа как маски. На языке битовых операций это условие можно записать как i | j = i, иными словами, $j$ — это число $i$, в котором занулили некоторые единичные биты.
Вычисление массива сумм по подмножествам является составным блоком для решения многих задач. К примеру, если элементы массива $a$ — это нули и единицы, обозначающие, является ли данное множество выигрышным, то элемент массива $b$, соответствующий какому-то множеству, будет положительным тогда и только тогда, когда какое-то его подмножество является выигрышным. Более того, значение в массиве $b$ — это количество выигрышных подмножеств данного множества.
Очевидный алгоритм для поиска массива $b$, работающий за время $O(n^2)$, — для каждой маски $i$ перебрать все другие маски и просуммировать те, которые являются подмасками $i$. Если воспользоваться алгоритмом перебора всех подмасок данной маски вместо того, чтобы перебирать все остальные маски в принципе, можно достичь асимптотики $O(3^k) = O(n^{\log_2 3}) \approx O(n^{1.585})$. Однако далее при помощи динамического программирования мы достигнем асимптотики $O(n \log n)$.
Давайте введем динамику dp, такую что dp[mask][ind] — это сумма по тем подмножествам маски $mask$, в которых изменялись только первые $ind$ битов. Если $ind = 0$, то мы не меняли никаких битов, и это просто $a_{mask}$. Далее при переходе от $ind$ к $ind + 1$ нужно, возможно, поменять текущий бит, либо не менять его. То есть, если бит $ind$ в $mask$ равен нулю, то мы ничего поменять не можем:
dp[mask][ind + 1] = dp[mask][ind]
Если же он равен единице, то мы его можем занулить:
dp[mask][ind + 1] = dp[mask][ind] + dp[mask ^ (1 << ind)][ind]
В конце dp[mask][k] как раз таки будет равно $b[mask]$, потому что это сумма по подмножествам, где мы можем занулять любые биты. Данную динамику можно слегка упростить, избавившись от второй размерности, и пересчитывать dp[mask][ind] в dp[mask][ind + 1] in-place. Итоговый алгоритм будет выглядеть следующим образом:
vector<int> sum_over_subsets(const vector<int>& a) {
vector<int> dp = a;
for (size_t bit = 1; bit < a.size(); bit <<= 1) {
for (size_t mask = 0; mask < a.size(); mask++) {
if ((mask & bit) != 0) {
dp[mask] += dp[mask ^ bit];
}
}
}
return dp;
}
Здесь переменная bit равна 1 << ind.
Очевидно, данный код работает за $O(n \cdot k) = O(n \log n)$, потому что первый цикл пробегает по всем $k$ битам, а второй — по всем $n$ маскам.
Обращение суммы по подмножествам
Иногда в задаче бывает нужно делать обратную операцию: по массиву сумм подмножеств $b$ построить изначальный массив $a$. Есть два подхода к придумыванию решения этой задачи. Можно придумать аналогичную динамику, в которой мы постепенно фиксируем биты. Но можно поступить проще. Можно воспользоваться общей идеей обращения алгоритмов, которая применима в разных контекстах: если у нас уже есть алгоритм, который делает какое-то прямое преобразование, то, чтобы сделать обратное преобразование, достаточно представить, что мы обращаем время вспять, и обратить все действия, которые мы сделали, в обратном порядке. Тогда нетрудно понять, что мы получим следующий алгоритм:
vector<int> inverse_sum_over_subsets(const vector<int>& b) {
vector<int> dp = b;
for (size_t bit = b.size() >> 1; bit >= 1; bit >>= 1) {
for (size_t mask = b.size() - 1; mask >= 0; mask--) {
if ((mask & bit) != 0) {
dp[mask] -= dp[mask ^ bit];
}
}
}
return dp;
}
Мы обратили порядок обоих циклов, тем самым выполняя все операции в обратном порядке, и если раньше мы что-то прибавили, то, чтобы обратить это действие, достаточно вычесть. Нетрудно заметить, что операции внутреннего цикла на одной итерации внешнего друг от друга не зависят, так что можно перевернуть внутренний цикл обратно. Аналогично, перебор битов во внешнем цикле по возрастанию был произвольным выбором, и для достижения того же результата мы можем перебирать биты в любом порядке, так что на самом деле оба цикла можно перевернуть обратно, и единственным отличием от прямого преобразования будет то, что вместо прибавления мы вычитаем:
vector<int> inverse_sum_over_subsets_simplified(const vector<int>& b) {
vector<int> dp = b;
for (size_t bit = 1; bit < b.size(); bit <<= 1) {
for (size_t mask = 0; mask < b.size(); mask++) {
if ((mask & bit) != 0) {
dp[mask] -= dp[mask ^ bit];
}
}
}
return dp;
}
И теперь уже совсем не составляет труда заметить, что получившийся код, — это в точности то, что мы бы получили, если бы попытались построить аналогичную динамику с последовательным фиксированием битов. Однако идея обращения всех операций в обратном порядке — это максимально общая концепция, которая применима во многих других контекстах.
Сумма по надмножествам
А что, если нам хочется найти не сумму по подмножествам, а сумму по надмножествам? Одно можно преобразовать в другое, заменив все нули в масках на единицы, а единицы на нули. Тогда зануление нескольких единиц в новой маске соответствует превращению нескольких нулей в единицы в старой маске, то есть переходу к какой-то надмаске. То есть нам надо просто поменять в алгоритме то, как мы обращаемся с нулями и с единицами. Если раньше, когда в маске какой-то бит был единицей, мы могли поменять его на ноль, то теперь наоборот: если какой-то бит в маске равен нулю, то мы можем поменять его на единицу. Тем самым, единственное, что нам надо поменять в коде, — это неравенство нулю на равенство:
vector<int> sum_over_supersets(const vector<int>& a) {
vector<int> dp = a;
for (size_t bit = 1; bit < a.size(); bit <<= 1) {
for (size_t mask = 0; mask < a.size(); mask++) {
if ((mask & bit) == 0) {
dp[mask] += dp[mask ^ bit];
}
}
}
return dp;
}
И обратная операция по получению массива $a$ по массиву $b$ аналогично получается заменой плюса на минус:
vector<int> inverse_sum_over_supersets(const vector<int>& b) {
vector<int> dp = b;
for (size_t bit = 1; bit < b.size(); bit <<= 1) {
for (size_t mask = 0; mask < b.size(); mask++) {
if ((mask & bit) == 0) {
dp[mask] -= dp[mask ^ bit];
}
}
}
return dp;
}
Сумма по подмножествам для проверки леммы Холла
Лемма Холла — стандартный способ проверки того, что в двудольном графе существует паросочетание, насыщающее левую долю. Она утверждает, что паросочетание есть тогда и только тогда, когда любое подмножество вершин левой доли соединено с не меньшим количество вершин правой доли. В общем случае, если размер левой доли равен $n$, для проверки потребуется минимум $O(2^n)$ времени, чтобы перебрать все подмножества левой доли. Легко заметить, к примеру, что нас интересуют только в некотором смысле максимальные по включению подмножества вершин левой доли. Иными словами, если к подмножеству вершин левой доли можно добавить еще одну вершину, не увеличив множество их соседей в правой доле, то если лемма Холла не выполнялась, то она и не начнет выполняться для нового подмножества. Таким образом, вместо подмножеств левой доли можно перебирать подмножества правой доли и смотреть на все вершины левой, которые соединены только с какими-то вершинами из этого подмножества. В общем случае это будет работать примерно за $O(2^m)$, где $m$ — количество вершин в правой доле, что совершенно не лучше, чем $O(2^n)$, ведь если $m < n$, паросочетания точно не существует. Однако в конкретных случаях можно сделать проверку сильно быстрее. К примеру, если у каждой вершины графа есть вес $w_v$, который означает что вершины $v$ есть $w_v$ копий. В таком случае реальное количество вершин в каждой доле равно сумме весов, однако перебирать нам нужно только подмножества типов вершин. В таком случае $m$ может быть сильно меньше $n$, и нам остается перебрать все подмножества правой доли и для каждого из них понять, сколько вершин левой доли с ними соединены.
Наивный алгоритм сделает это за количество ребер в графе для каждого из $2^m$ подмножеств, чтобы найти все вершины левой доли, которые соединены только с этим подмножеством, однако при помощи сумм по подмножествам можно сделать это сильно быстрее. Определим a[mask] как сумму весов всех вершин левой доли, которые соединены ровно с вершинами правой доли из маски. Такой массив можно посчитать за линейное время от количества ребер в графе. Если мы теперь за $O(2^m \cdot m)$ построим массив b[mask], являющийся суммой по подмножествам массива a, то в нем как раз таки и будут храниться суммарные веса всех вершин левой доли, которые соединены с каким-то подмножеством вершин из маски. Теперь лишь остается для каждой маски сравнить это число с суммой весов вершин из маски, и тем самым мы проверим лемму Холла. К примеру, такое решение можно использовать, когда $n$ имеет порядок $10^5$, а $m \approx 20$.
or-свертка и and-свертка
Свертки последовательностей — часто встречающийся шаг в решении задач. К примеру, быстрое преобразование Фурье (знание быстрого преобразования Фурье необязательно для понимания дальнейших алгоритмов) позволяет по двум массивам $a$ и $b$ считать массив $c$, такой что $c_i = \sum_{j=0}^i a_j \cdot b_{i-j}$ за $O(n \log n)$, хотя тривиальный алгоритм сделал бы это только за $O(n^2)$. Такая свертка двух последовательностей бывает очень полезной во многих задачах. Альтернативно можно записать ее так:
$$c_i = \sum_{j, k : j + k = i} a_j \cdot b_{k}$$
То есть мы суммируем произведения элементов с индексами, в сумме дающими $i$. Однако по аналогии можно рассмотреть и другую свертку, которая тоже бывает полезна при решении задач: or-свертку. Формула аналогична:
$$c_i = \sum_{j, k : j | k = i} a_j \cdot b_k$$
Найти последовательность $c$ по последовательностям $a$ и $b$ за $O(n \log n)$ нам как раз поможет сумма по подмножествам. Для простоты увеличим длины массивов $a$ и $b$ до ближайшей степени двойки, и тогда мы можем считать, что все три массива имеют длину $n = 2^k$. По аналогии с быстрым преобразованием Фурье, мы сначала преобразуем последовательности $a$ и $b$ в альтернативный вид, в котором для получения свертки достаточно перемножить их поэлементно, а затем сделаем обратное преобразование и тем самым получим последовательность $c$. Для быстрого преобразования Фурье это преобразование было вычислением значений многочленов в фиксированных точках, а в нашем случае это будет переход от последовательностей к их суммам по подмножествам. За $O(2^k \cdot k) = O(n \log n)$ мы можем построить последовательности $f$ и $g$, которые являются суммами по подмножествам $a$ и $b$ соответственно. Тогда утверждается, что последовательность $h$, полученная поэлементным произведением $f$ и $g$ ($h_i = f_i \cdot g_i$) является суммой по подмножествам последовательности $c$. Почему это так? По определению суммы по подмножествам, мы можем расписать $h_i$ следующим образом:
$$h_i = f_i \cdot g_i = \sum_{j : j \subset i} a_j \cdot \sum_{k : k \subset i} b_k$$
Произведение сумм — это сумма попарных произведений, так что
$$h_i = \sum_{j : j \subset i, k : k \subset i} a_j \cdot b_k$$
То, что $j$ и $k$ являются подмасками маски $i$ эквивалентно тому, что их побитовый or является подмаской $i$, поэтому эту сумму можно переписать следующим образом:
$$h_i = \sum_{j, k : j | k \subset i} a_j \cdot b_k$$
А это в точности и есть сумма по подмножествам последовательности $c$, ведь в последовательности $c$ нас интересует сумма произведений по всем парам маскок, у которых or равен $i$, а здесь мы берем сумму произведений по всем парам масок, у которых or — это подмаска $i$. Иными словами, если мы обозначим $j | k$ за $l$, мы можем переписать сумму выше как
$$h_i = \sum_{l \subset i} \sum_{j, k : j | k = l} a_j \cdot b_k = \sum_{l \subset i} c_l$$
Таким образом, чтобы из последовательности $h$ получить последовательность $c$, достаточно применить алгоритм обращения суммы по подмножествам. Итоговый алгоритм or-свертки выглядит следующим образом:
// a.size() == b.size() == 2^k
vector<int> or_convolution(const vector<int>& a, const vector<int>& b) {
vector<int> f = sum_over_subsets(a);
vector<int> g = sum_over_subsets(b);
vector<int> h(f.size());
for (size_t i = 0; i < f.size(); i++) {
h[i] = f[i] * g[i];
}
vector<int> c = inverse_sum_over_subsets(h);
return c;
}
and-свертка
Легко понять, что для операции and все работает точно так же. Если мы хотим получить свертку
$$c_i = \sum_{j, k : j \& k = i} a_j \cdot b_k$$
достаточно взять сумму по надмножествам последовательностей $a$ и $b$, перемножить результаты, а потом взять обращение суммы по надмножествам для этого произведения. Для доказательства можно либо провести абсолютно аналогичные рассуждения, либо вспомнить, что сумма по надмножествам отличалась от суммы по подмножествам заменой в уме всех битов $1$ на $0$ и наоборот. Операции or и and отличаются друг от друга ровно такой же заменой.
xor-свертка
Преобразование Уолша — Адамара и его применение
Мы разобрались с or-сверткой и and-сверткой. Остается разобраться с последней битовой операцией — xor. Иными словами, по двум последовательностям $a$ и $b$ мы хотим построить последовательность $c$, такую что
$$c_i = \sum_{j, k : j \oplus k = i} a_j \cdot b_k = \sum_j a_j \cdot b_{i \oplus j}$$
где $\oplus$ — это операция xor. Для этого нам опять нужно придумать такое преобразование, чтобы поэлементно перемноженные преобразования от $a$ и $b$ дали преобразование от $c$. Мы определим его следующим образом. Преобразованием последовательности $d$ назовем последовательность $e$, такую что
$$e_i = \sum_{j} (-1)^{|i \& j|} d_j$$
Иными словами, это знакочередующаяся сумма всех элементов $d$, где маска $j$ берется с плюсом, если их пересечение с маской $i$ четного размера, и с минусом, если нечетного. Такое преобразование называется преобразованием Уолша — Адамара. Почему же такое преобразование подходит? Это следует из того, что четность суммы размеров двух масок равна четности размера их ксора (симметрической разности). Обозначим за $f$ преобразование последовательности $a$, за $g$ преобразование последовательности $b$ и за $h$ поэлементное произведение $f$ и $g$ ($h_i = f_i \cdot g_i$). Мы хотим доказать, что $h$ — это преобразование Уолша — Адамара от $c$. Распишем:
$$h_i = f_i \cdot g_i = \sum_{j} (-1)^{|i \& j|} a_j \cdot \sum_{k} (-1)^{|i \& k|} b_k$$
Как и в случае с or-сверткой, распишем произведение сумм как сумму попарных произведений:
$$h_i = \sum_{j, k} (-1)^{|i \& j| + |i \& k|} a_j \cdot b_k$$
Число $(-1)^{|i \& j| + |i \& k|}$ зависит от четности $|i \& j| + |i \& k|$. На картинке ниже эта сумма — размер синей области плюс размер зеленой области плюс дважды размер голубой. Но дважды учтенная голубая область не влияет на четность, поэтому четность суммы $|i \& j| + |i \& k|$ равна четности $|i \& (j \oplus k)|$ (размер пересечения $i$ и симметрической разности $j$ и $k$), потому что это в точности сумма размеров синей и зеленой областей, что легко понять по картинке.
Таким образом, значение $h_i$ можно переписать:
$$h_i = \sum_{j, k} (-1)^{|i \& (j \oplus k)|} a_j \cdot b_k$$
Обозначим $l := j \oplus k$ и просуммируем отдельно по каждому значению $l$:
$$h_i = \sum_{l} \sum_{j, k : j \oplus k = l} (-1)^{|i \& l|} a_j \cdot b_k$$
Член $(-1)^{|i \& l|}$ можно вынести из внутренней суммы, после чего внутренняя сумма в точности будет равна $c_l$ по определению:
$$h_i = \sum_l (-1)^{|i \& l|} \sum_{j, k : j \oplus k = l} a_j \cdot b_k = \sum_l (-1)^{|i \& l|} c_l$$
Таким образом, мы получили, что $h_i$ в точности равно преобразованию Уолша — Адамара для последовательности $c$. Поэтому, чтобы найти последовательность $c$, нам остается лишь понять, как эффективно производить это преобразование и его обращение.
Вычисление преобразования Уолша — Адамара и его обращения
Мы хотим посчитать по последовательности $a$ такую последовательность $f$, что
$$f_i = \sum_{j} (-1)^{|i \& j|} a_j$$
По аналогии с суммой по подмножествам, мы будем считать последовательность $f$ при помощи динамики, постепенно перебирая биты, разрешая им меняться. Введем динамику dp, такую что dp[mask][ind] — это знакочередующаяся сумма по тем подмножествам маски $mask$, в которых изменялись только первые $ind$ битов, и знак определяется размером пересечения масок среди первых $ind$ битов. Если $ind = 0$, то мы не меняли никаких битов, и это просто $a_{mask}$. Далее, при переходе от $ind$ к $ind + 1$ нужно, возможно, поменять текущий бит, либо не менять его. Если текущий бит в маске равен нулю, то меняем мы его или нет, знаки остаются теми же самыми, потому что размер пересечения масок не изменился:
dp[mask][ind + 1] = dp[mask][ind] + dp[mask ^ (1 << ind)][ind]
Если же текущий бит равен единице, и мы его не меняем, то в пересечение добавляется один новый бит, и знак меняется:
dp[mask][ind + 1] = dp[mask ^ (1 << ind)][ind] - dp[mask][ind]
В конце dp[mask][k] как раз таки будет равно $f[mask]$, потому что это знакочередующаяся сумма по подмножествам, где мы можем менять любые биты. Данную динамику можно слегка упростить, избавившись от второй размерности, и пересчитывать dp[mask][ind] в dp[mask][ind + 1] in-place. Однако в данном случае нужно быть осторожным. Два значения, отличающиеся на (1 << ind), в данном случае оба меняются одновременно, так что необходимо обновлять их одномоментно, чтобы значения, через которые мы пересчитываемся, действительно были равны dp[...][ind] Итоговый алгоритм будет выглядеть следующим образом:
vector<int> hadamard_transform(const vector<int>& a) {
vector<int> dp = a;
for (size_t bit = 1; bit < a.size(); bit <<= 1) {
for (size_t mask = 0; mask < a.size(); mask++) {
if ((mask & bit) == 0) {
int u = dp[mask], v = dp[mask ^ bit];
dp[mask] = u + v;
dp[mask ^ bit] = u - v;
}
}
}
return dp;
}
Здесь переменная bit равна 1 << ind.
Очевидно, что данный код работает за $O(n \cdot k) = O(n \log n)$, потому что первый цикл пробегает по всем $k$ битам, а второй — по всем $n$ маскам.
Чтобы получить обратное преобразование, можно либо посмотреть на формулу определения преобразования и понять, что если мы применим его к $f$, то мы получим изначальную последовательность $a$, в которой все элементы умножены на $n$, либо же воспользоваться общей стратегией: обратить все действия в обратном порядке. В этом случае это немного сложнее, чем раньше. Перед нами есть значения $x = u + v$ и $y = u - v$, и мы хотим обратно получить $u$ и $v$. Несложно понять, что $u = (x + y) / 2$ и $v = (x - y) / 2$. Также нетрудно понять, что циклы (по аналогии с суммой по подмножествам) переворачивать необязательно:
vector<int> inverse_hadamard_transform(const vector<int>& f) {
vector<int> dp = f;
for (size_t bit = 1; bit < f.size(); bit <<= 1) {
for (size_t mask = 0; mask < f.size(); mask++) {
if ((mask & bit) == 0) {
int x = dp[mask], y = dp[mask ^ bit];
dp[mask] = (x + y) / 2;
dp[mask ^ bit] = (x - y) / 2;
}
}
}
return dp;
}
Обратите внимание, что отличие между прямой и обратной сверткой заключается только в том, что в обратной мы делим пополам. Можно заметить, что мы поделим пополам логарифм раз, так что на самом деле, чтобы сделать обратное преобразование, достаточно сделать прямое, а в конце поделить все элементы на размер массива.
Главное только помнить, что во всех этих преобразованиях можно менять порядки прохода по циклам, но не циклы местами!
Наконец, итоговая xor-свертка будет выглядеть абсолютно аналогично or-свертке:
// a.size() == b.size() == 2^k
vector<int> xor_convolution(const vector<int>& a, const vector<int>& b) {
vector<int> f = hadamard_transform(a);
vector<int> g = hadamard_transform(b);
vector<int> h(f.size());
for (size_t i = 0; i < f.size(); i++) {
h[i] = f[i] * g[i];
}
vector<int> c = inverse_hadamard_transform(h);
return c;
}
Альтернативный математический взгляд
Если вы знакомы с алгеброй, на xor-свертку можно смотреть иначе. Преобразование Фурье воспринимает входную последовательность как коэффициенты многочлена от одной переменной и считает значения этого многочлена в каких-то точках. Преобразование Уолша — Адамара же воспринимает входную последовательность как коэффициенты многочлена от $k$ переменных, где в каждый одночлен каждая переменная входит в степени $0$ или $1$ (единичные биты маски в точности и говорят, какие переменные стоят перед этим коэффициентом). К примеру, последовательность $a_0, a_1, a_2, a_3$ воспринимается как многочлен $P(x_1, x_2) = a_0 + a_1 x_1 + a_2 x_2 + a_3 x_1 x_2$. Преобразование Уолша — Адамара в точности находит значения этого многочлена во всех $2^k$ точках, где все переменные равны $\pm 1$. А так как $x_i^2 = 1$ в таком случае, то произведение двух многочленов $P$ и $Q$ от последовательностей $a$ и $b$ как раз таки дает многочлен $R$ для свертки этих последовательностей $c$.
Сумму по подмножествам можно альтернативно воспринимать таким же образом. Только считаем мы значения в точках, где значения всех переменных — это $0$ или $1$. В таком случае выполняется условие $x_i^2 = x_i$. Для and-свертки, как всегда, достаточно поменять нули и единицы местами.
Источники
- Пункт 15 этого блога рассказывает об альтернативном математическом взгляде на свертки.
- Старая более алгебраическая версия этой главы.