Решаем систему логических уравнений на Rust
• Павел Никитин • руководства • поддержите на Patreon
Язык программирования Rust сосредоточен на безопасности, скорости и параллелизме. Давайте проверим, насколько Rust быстрый в сравнении с другими компилируемыми языками программирования C и Vala в решении системы логических уравнений. Для этого реализуем алгоритм решения интересной головоломки Lights Off на языках
- Rust
- C
- Vala
и сравним их скорость. Подробно рассмотрим реализацию на языке Rust, реализации на других языках можно изучить в открытом репозитории.
Задача
Дано прямоугольное поле, каждая ячейка которого включена или выключена. При переключении ячейки переключаются также её недиагональные соседи. Требуется выключить все ячейки за минимальное количество переключений.
Решение
По заданному полю формируем систему логических уравнений с числом неизвестных, совпадающим с числом ячеек. Решаем систему методом Гаусса. Если ранг системы равен числу неизвестных, то получаем единственное решение. Иначе проверяем систему на противоречие. Если противоречия нет, то выбираем решение с минимальным весом среди подмножеств остатка за рангом системы.
Реализация
Обзор
Загружаем исходное поле в битовую матрицу. Битовая матрица является вектором битовых векторов. Битовый вектор хранится в векторе целых чисел без знака, разрядность которых совпадает с разрядностью платформы. По битовой матрице поля формируем битовую матрицу системы уравнений. Решаем систему методом Гаусса, выбираем решение с минимальным весом. Преобразуем битовый вектор решения в битовую матрицу результата и выводим её.
Битовый вектор
Объявляем структуру, владеющую вектором целых чисел без знака usize, разрядность которых совпадает с разрядностью платформы. Такой подход позволяет компактно хранить битовый вектор и оперировать множеством битов.
1 2 3 4 | pub struct BitVec { buf: Vec<usize>, len: usize, } |
Реализуем чтение и запись отдельных битов вектора. Компактные методы объявляем с директивой #[inline] для избежания вызовов и расширения манёвра для оптимизации использования регистров процессора в результирующем машинном коде.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | #[inline] pub fn get(&self, index: usize) -> bool { (self.buf[buf_index(index)] >> bit_index(index)) & 1 == 1 } #[inline] pub fn set(&mut self, index: usize, value: bool) { if value { self.buf[buf_index(index)] |= bit_mask(index); } else { self.buf[buf_index(index)] &= !bit_mask(index); } } |
Вспомогательные функции для битового доступа должны быть максимально быстрыми, поскольку они работают в циклах с большим количеством итераций.
1 2 3 4 5 6 7 8 9 | #[inline] fn bit_index(index: usize) -> usize { index & BIT_MASK } #[inline] fn bit_mask(index: usize) -> usize { 1usize << bit_index(index) } |
Битовая маска для доступа к отдельным битам зависит от разрядности платформы, которую можно получить из параметра target_pointer_width макроса cfg.
1 2 3 4 5 6 7 8 9 10 11 | #[cfg(target_pointer_width = "32")] pub const WORDSIZE: usize = 32; #[cfg(target_pointer_width = "32")] const BIT_SHIFT: u8 = 5; #[cfg(target_pointer_width = "64")] pub const WORDSIZE: usize = 64; #[cfg(target_pointer_width = "64")] const BIT_SHIFT: u8 = 6; const BIT_MASK: usize = WORDSIZE-1; |
Исходный код битового вектора: bitvec.rs. Заметим, что добавление тестов в комментариях к методам существенно облегчает жизнь при тестировании и позволяет документировать их.
Битовая матрица
Структура битовой матрицы владеет вектором битовых векторов.
1 2 3 | pub struct BitMat { rows: Vec<BitVec> } |
В целом её реализация не нуждается в подробных комментариях, за исключением двух методов. Это метод swap для перестановки двух строк матрицы и xor для побитового сложения двух строк. Сложность в том, что эти методы одновременно считывают и изменяют состояние битовой матрицы, а это противоречит концепции владения ресурсами в Rust. Чтобы успокоить компилятор, заключим код этих методов в блок unsafe и используем сырые указатели.
1 2 3 4 5 6 7 8 | #[inline] pub fn swap(&mut self, row_i: usize, row_j: usize) { unsafe { let p_row_i = &mut self.rows[row_i] as *mut BitVec; let p_row_j = &mut self.rows[row_j] as *mut BitVec; ptr::swap(p_row_i, p_row_j); } } |
В методе xor мы получаем дивиденды от затеи с битовыми векторами. Многократно повторяющаяся операция сложения двух строк в методе Гаусса, адаптированном для поля по модулю два, происходит быстро за счёт сокращения количества операций, кратного разрядности платформы. Так, на 64-битной платформе одной инструкцией процессора мы складываем 64 элемента системы уравнений!
1 2 3 4 5 6 7 8 9 10 11 | #[inline] pub fn xor(&mut self, row_i: usize, row_j: usize) { unsafe { let p_row_j = self.rows[row_j].buf() as *const Vec<usize>; let len = self.rows[row_j].buf().len(); for i in 0..len { self.rows[row_i].buf_mut()[i] ^= (&*p_row_j)[i]; } } } |
Исходный код битовой матрицы: bitmat.rs.
Решение системы логических уравнений
Применим наш битовый опыт в решении системы логических уравнений методом Гаусса. Создадим структуру алгоритма решения системы, владеющего системой, её рангом, числом решений и минимальным весом решения.
1 2 3 4 5 6 | pub struct BitGauss { sys: BitMat, rank: usize, n_solutions: usize, min_weight: usize, } |
Метод Гаусса позволяет преобразовать матрицу системы к единичной матрице, содержащей единицы на главной диагонали.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 | pub fn gauss(&mut self) { let n_rows = self.sys.n_rows(); // Преобразуем левую часть матрицы системы к единичной матрице for i in 0..n_rows { // Найдём единицу в столбце `i` и переставим местами найденную строку со строкой `i` for j in i..n_rows { if self.sys.get(j, i) { self.sys.swap(j, i); } } // Пропускаем столбец, если он не содержит единицу if !self.sys.get(i, i) { continue; } // Увеличиваем ранг системы; столбец содержит единицу self.rank = i+1; // Обнуляем элементы столбца `i` кроме строки `i` for j in 0..n_rows { if self.sys.get(j, i) && j != i { self.sys.xor(j, i); } } } } |
Рассмотрим алгоритм выбора решения с минимальным весом после применения метода Гаусса к системе уравнений. Если ранг системы равен числу неизвестных, то получаем единственное решение. Иначе проверяем систему на противоречие. Если противоречия нет, то выбираем решение с минимальным весом среди подмножеств остатка за рангом системы. Для перебора подмножеств остатка запустим индекс решения от нуля до двух в степени разницы между числом неизвестных и рангом системы. Двоичное представление значения индекса решения в битовом векторе даёт текущее подмножество остатка: единица — берём элемент, ноль — не берём элемент в подмножество. Для каждой строки до ранга системы включительно аккумулируем элементы из столбцов подмножества остатка с элементом правой части системы. Подсчитываем вес решения как количество единиц в аккумуляторе и подмножестве остатка. Если нашли решение с меньшим весом, то собираем вектор решения из аккумулятора и подмножества остатка. Это ключевой алгоритм в решении головоломки, приведём его исходный код.
pub fn solve(&mut self) -> Option<BitVec> {
// Применяем метод Гаусса к системе уравнений
self.gauss();
// Проверяем систему на противоречие
let n_rows = self.sys.n_rows();
let n_cols = self.sys.n_cols();
let n_vars = n_cols-1;
let rank = self.rank;
// Если есть противоречие, то система не имеет решения
for i in rank..n_rows {
if self.sys.get(i, n_vars) {
return None
}
}
// Создаём битовый вектор решения
let mut solution = BitVec::with_length(n_vars);
// Если ранг системы совпадает с числом неизвестных,
// то собираем решение из правой части системы
if rank == n_vars {
for i in 0..n_rows {
solution.set(i, self.sys.get(i, n_vars));
}
}
// Иначе выбираем решение с минимальным весом из множества остатка за рангом системы
else {
let n_rest = n_vars - rank;
self.n_solutions = 1usize << n_rest;
self.min_weight = n_cols;
let mut weight;
let mut rest = BitVec::with_length(n_rest);
let mut accumulator = BitVec::with_length(n_rows);
// Ищем решение с минимальным числом единиц
for i in 0..self.n_solutions {
// Создаём подмножество остатка по двоичному представлению индекса решения
rest.buf_mut()[0] = i;
// Обнуляем аккумулятор
accumulator.setall(false);
// Аккумулируем решение с текущим индексом
// из подмножества остатка и правой части системы
for j in 0..rank {
for k in 0..n_rest {
if rest.get(k) {
accumulator.xor(j, self.sys.get(j, rank + k));
}
}
accumulator.xor(j, self.sys.get(j, n_vars));
}
// Взвешиваем аккумулятор и подмножество остатка
weight = (accumulator.count_ones() + rest.count_ones()) as usize;
// Конструируем решение с меньшим весом
if weight < self.min_weight {
self.min_weight = weight as usize;
// Берём первую часть элементов вектора решения из аккумулятора
for j in 0..rank {
solution.set(j, accumulator.get(j));
}
// Берём вторую часть элементов вектора решения из подмножества остатка
for j in 0..n_rest {
solution.set(rank + j, rest.get(j));
}
}
}
}
Some(solution)
}
Исходный код алгоритма решения системы уравнений: bitalg.rs.
Построение системы уравнений
Применим наш опыт решения систем логических уравнений для решения головоломки. Создадим структуру решателя головоломки, владеющего алгоритмом решения системы и размером поля.
pub struct LightsSolver {
alg: BitGauss,
n_rows: usize,
n_cols: usize,
}
В конструкторе решателя опишем как построить систему уравнений по заданному полю головоломки. Поскольку у нас каждая ячейка — неизвестная переменная в системе уравнений, то число уравнений определяется произведением числа строк и столбцов поля. Число элементов системы является квадратом числа ячеек вместе с вектором правой части. Здесь мы получаем ещё один гешефт от затеи с битовыми векторами. Для поля 100 на 100 ячеек количество элементов системы будет 100 010 000! Логический тип обычно имеет размер разрядности платформы. Так, для 64-битной платформы наша система весила бы 763 Мб вместо 12 Мб! Вернёмся к построению системы. Владение полем решателю не требуется, поэтому мы на время построения заимствуем ссылку. Помня о том, что при переключении ячейки переключаются её недиагональные соседи, обойдём каждую ячейку поля и её соседей справа, слева, сверху и снизу, определяя индекс ячейки в матрице системы уравнений.
#[inline]
fn sys_index(row: isize, col: isize, n_rows: isize, n_cols: isize) -> isize {
if 0 <= row && row < n_rows as isize && 0 <= col && col < n_cols as isize {
n_cols * row + col
} else {
-1
}
}
Если индекс ячейки положительный, то устанавливаем соответствующий элемент в левой части матрицы системы. Правую часть матрицы устанавливаем по полю головоломки в зависимости от того, включена там ячейка или нет.
pub fn from(field: &BitMat) -> LightsSolver {
let n_rows = field.n_rows() as isize;
let n_cols = field.n_cols() as isize;
let n = n_rows * n_cols;
let mut himself: isize;
let mut neighbor: isize;
let mut sys = BitMat::with_size(n as usize, (n + 1) as usize);
for row in 0..n_rows {
for col in 0..n_cols {
himself = sys_index(row, col, n_rows, n_cols);
sys.set(himself as usize, himself as usize, true);
neighbor = sys_index(row, col + 1, n_rows, n_cols);
if neighbor >= 0 {
sys.set(himself as usize, neighbor as usize, true);
}
neighbor = sys_index(row, col - 1, n_rows, n_cols);
if neighbor >= 0 {
sys.set(himself as usize, neighbor as usize, true);
}
neighbor = sys_index(row + 1, col, n_rows, n_cols);
if neighbor >= 0 {
sys.set(neighbor as usize, himself as usize, true);
}
neighbor = sys_index(row - 1, col, n_rows, n_cols);
if neighbor >= 0 {
sys.set(neighbor as usize, himself as usize, true);
}
sys.set(himself as usize, n as usize, field.get(row as usize, col as usize));
}
}
let alg = BitGauss::from(sys);
LightsSolver {
alg: alg,
n_rows: n_rows as usize,
n_cols: n_cols as usize,
}
}
Почти всё готово. Нам осталось решить полученную систему и преобразовать вектор решения в матрицу ответов, в которой единицы показывают, какие ячейки нужно переключать в исходном поле в произвольном порядке, чтобы полностью его выключить. Lights Off.
pub fn solve(&mut self) -> Option<BitMat> {
match self.alg.solve() {
Some(syssol) => {
let mut sol = BitMat::with_size(self.n_rows, self.n_cols);
for row in 0..self.n_rows {
for col in 0..self.n_cols {
if syssol.get(self.n_cols * row + col) {
sol.set(row, col, true);
}
}
}
Some(sol)
},
None => None
}
}
Исходный код решателя головоломки: ligsol.rs. В файле main.rs реализован незатейливый консольный интерфейс.
Производительность
Условия тестирования
Выберем несколько размеров квадратного поля с изначально включёнными ячейками. Здесь нужно обратить внимание на количество решений для заданного размера поля. Для поля 30×30 нужно перебрать более миллиона решений, для поля 50×50 существует 256 решений, для полей 100×100, 150×150, 200×200 существует по одному решению. Тестировать будем на платформе со следующими параметрами:
- Процессор: Intel Core i5-7200U
- Максимальная частота: 3.10 ГГц
- Кэш процессора: 3 Мб
- Оперативная память: 16 Гб
- Операционная система: Arch Linux
- Ядро: 4.10.13-1-ARCH
- Toolchain: nightly-x86_64-unknown-linux-gnu
Результат
Сведём результаты тестирования в таблицу, отражающую время решения задачи заданного размера в секундах для разных языков программирования. Флаг сборки для Rust установили --release, для C -O3, для Vala -O3.
| Поле | Rust | C | Vala |
|---|---|---|---|
| 50×50 | 0 | 0 | 3 |
| 100×100 | 4 | 4 | 166 |
| 150×150 | 50 | 45 | 1835 |
| 200×200 | 215 | 210 | 9745 |
Анализ
Алгоритм реализован на трёх компилируемых языках программирования, скомпилированные программы работают на одном и том же процессоре и операционной системе в условиях одинаковой нагрузки. Реализация на языке Vala оказалась существенно медленней реализаций на языках C и Rust. Анализ промежуточного исходного кода на языке C для реализации на языке Vala показал отсутствие выделения и освобождения памяти в циклах с большим количеством итераций, но обнаружил громоздкий код с большим количеством лишних инструкций, в результате чего результирующий машинный код оказался неэффективным. Реализация на языке Rust получилась чуть медленнее реализации на языке C, что говорит нам о высоком качестве результирующего машинного кода компилятора Rust, сопоставимого с машинным кодом компилятора C.
Перспективы
Метод Гаусса плохо поддаётся распараллеливанию, но можно хорошо распараллелить перебор решений. Это делает язык Rust приоритетным кандидатом из-за его выразительных возможностей распараллеливания. Можно прикрутить графический интерфейс и получится забавная игрушка.
Заключение
Rust оправдывает звание скоростного языка программирования!
