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

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

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

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

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

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

Работать с ним нужно вот так:

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

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

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

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

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

Обычно его решают так:

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

…а с битсетом оно разгоняется так:

bitset<lim> 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\), используя ровно \(n\) переходов. Иногда нам не нужно знать число способов — нам хватит знания, можно ли вообще через \(n\) ходов там оказаться. Тогда вместо числового умножения нам хватит битового умножения, то есть &:

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

matrix operator*(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\) лампочек. Известно текущее состояние всех лампочек, нужно восстановить по нему состояние переключаетелей.

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

\[ \begin{cases} a_{11} x_1 + a_{12} x_2 + \ldots + a_{1n} x_n \equiv b_1 \pmod 2\\ a_{21} x_1 + a_{22} x_2 + \ldots + a_{2n} x_n \equiv b_2 \pmod 2\\ \ldots \\ a_{n1} x_1 + a_{n2} x_2 + \ldots + a_{nn} x_n \equiv b_n \pmod 2 \end{cases} \]

Здесь \(x\) — состояния переключателей, \(b\) — состояния лампочек, а матрица \(A\) содержит информацию о том, влияет ли переключатель на лампочку.

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

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\).