Битсет и битовое сжатие

Нумерованные множества можно представлять в виде булевых массивов: если элемент $x$ присутствует в множестве, то $x$-тый элемент массива будет равен единице, и нулю в противном случае.

Разные теоретико-множественные операции часто можно свести к поэлементным операциям булевыми массивами. Например:

• Объединение множеств — побитовое ИЛИ.
• Пересечение множеств — побитовое И.
• Симметрическая разность — побитовый XOR (исключающее ИЛИ).

Процессор устроен так, что работает не с отдельными битами, а сразу с блоками по, например, 32 или 64 бита — эта величина называется машинным словом — и поэтому операции, затрагивающие лишь один бит, на самом деле «стоят» столько же, сколько и операции над всеми битами int или long.

Здесь появляется следующая идея оптимизации: сгруппировать элементы булева массива в блоки размера 64 и каждый такой блок представить 64-битным двоичным числом. Тогда можно применять соответствующие побитовые операцию сразу к 64 элементам и тратить на это один процессорный такт вместо 64-х.

#std::bitset

Это всё несложно кодить и вручную, но в STL всё уже сделано до нас. bitset — структура, ведущая себя как большое двоичное число со всеми стандартными битовыми операциями:

const int lim = 1000;
bitset<lim> b; // создать битсет размера lim (должно быть константой)
b.set();       // заполнить единицами
b.reset();     // заполнить нулями
b.flip();      // заменить единички на нули и наоборот
b.count();     // посчитать число единичек
cout << b;     // вывести бинарную строку
cout << b[42]; // вывести 43-ий (нумерация с нуля, как у массивов) бит "справа"

Также для битсетов работает вся битовая арифметика: &, |, ^, ~, <<, >> и их варианты с [operator]=.

Примечание. Часто «асимптотику» использующих битовое сжатие алгоритмов пишут как $O(f / 64)$. Автору не нравится такая нотация, потому что численную константу формально можно сократить, и, вообще, на разных архитектурах и у разных реализаций будут разные константы. Вместо этого будем везде писать $O(f / w)$, где $w$ означает размер машинного слова.

#Задача о рюкзаке

Даны $n$ предметов с положительными целыми весами $a_i$ и рюкзак размера $m$. Требуется выбрать подмножество предметов с максимальной суммой, не превышающий размер рюкзака.

Обычное решение за $O(n \cdot m)$:

bool dp[m] = {}; // так можно его заполнить нулями
dp[0] = 1;
for (int i = 0; i < n; i++)
    for (int x = m - a[i]; x >= 0; x--)
        dp[x + a[i]] |= dp[x];

…битовым сжатием разгоняется до $O(n \cdot m / w)$ так:

bitset<m> b;
b[0] = 1;
for (int i = 0; i < n; i++)
    b |= b << a[i];

#Цикл длины 3

Нужно узнать, есть ли цикл длины 3 в ориентированном графе из $n$ вершин, заданном своей матрицей смежности.

Обернем матрицу смежности в массив битсетов, и тогда задача решается за $O(n^3 / w)$:

bitset<maxn> g[maxn]; // матрица смежности
for (int a = 0; a < n; a++) {
    for (int b = 0; b < n; b++) {
        if (g[a][b] && (~g[a] & g[b]).any()) {
            // цикл найден
        }
    }
}

Бенчмарк: на серверах CodeForces этот код при $n = 5000$ работает за 7 секунд.

#Перемножение матриц

В задачах на подсчет числа путей часто используется факт, что матрица смежности графа, возведенная в степень $n$, имеет комбинаторный смысл: на пересечении $a$-ой строки и $b$-того столбца матрицы $G^n$ будет записано количество способов дойти из вершины $a$ в вершину $b$, используя ровно $n$ переходов.

В некоторых задачах нам не нужно знать именно число способов дойти из $a$ в $b$ — нам достаточно знания, можно ли вообще там оказаться за $n$ ходов. Тогда вместо числового умножения нам хватит битового умножения, то есть &:

typedef bitset<maxn> t;
typedef array<t, maxn> matrix;

matrix matmul(matrix a, matrix b) {
    matrix c;
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            if (a[i][j])
                c[i] |= b[j];
    return c;
}

#Метод Гаусса

Иногда встречаются задачи, требующие решения системы линейных уравнений. Большую часть из них на самом деле можно решить над полем $\mathbb{Z}_2$ — «по модулю 2».

Есть $n$ переключателей лампочек. Каждый активированный переключатель меняет состояние (включает или выключает) какого-то подмножества из $n$ лампочек. Известно текущее состояние всех лампочек, требуется восстановить по нему состояние переключателей.

Нас по сути просят решить следующую систему:

Здесь $x$ — состояния переключателей, $b$ — состояния лампочек, а бит $a_{ij}$ отражает, влияет ли переключатель $i$ на лампочку $j$.

В таком случае можно значительно ускорить и упростить обычный метод Гаусса:

t gauss(matrix a) {
    for (int i = 0; i < n; i++) {
        int nonzero = i;
        for (int j = i+1; j < n; j++)
            if (a[j][i])
                nonzero = j;
        swap(a[nonzero], a[i]);
        for (int j = 0; j < n; j++)
            if (j != i && a[j][i])
                a[j] ^= a[i];
    }
    t x;
    for (int i = 0; i < n; i++)
        x[i] = a[i][n] ^ a[i][i];
    return x;
}

Код находит вектор $x$ из уравнения $Ax = b$ при условии, что решение существует и единственно. Для простоты кода предполагается, что вектор $b$ приписан справа к матрице $A$.

#Детали реализации

На самом деле, на высоких уровнях оптимизации (g++ -O3 -march=native ...) компилятор в конечном итоге будет использовать группы не по 64, а по 256 бит, используя SIMD-инструкции. Однако, для реализации .count() самое быстрое, что можно придумать — вызывать инструкцию popcnt по очереди от каждого 64-битного числа (в GCC она доступна как встроенный интринзик __builtin_popcount), поэтому этот метод работает в несколько раз медленнее побитовых операций.

Операции с битсетами довольно просто реализовывать и с нуля: нужен только цикл с достаточно простым телом, и компилятор сделает всё сам, никакой тёмной магии при этом не требуется.