Метод четырёх русских для умножения матриц

Постановка задачи

Рассмотрим следующую задачу: «Дано две квадратных матрицы [math]A_{[n \times n]}[/math] и [math]B_{[n \times n]}[/math], состоящие из нулей и единиц. Нужно найти их произведение. При этом, все операции выполняются по модулю [math]2[/math].»

Простое решение

Если мы будем считать произведение матриц [math]C = A \cdot B[/math] по определению([math]c_{i, j} = \sum\limits_{k = 1}^n a_{i,k}b_{k,j}[/math]), то трудоёмкость алгоритма составит [math]O(n^3)[/math] — каждый из [math]n^2[/math] элементов результирующей матрицы [math]C[/math] вычисляется за время, пропорциональное [math]n[/math].

Сейчас будет показано, как немного уменьшить это время.

Предподсчёт

Воспользуемся следующим приёмом. Возьмём некоторое целое число [math]k[/math]. Для всех возможных пар двоичных векторов длины [math]k[/math] подсчитаем и запомним их скалярное произведение по модулю [math]2[/math].

Сжатие матриц

Воспользуемся полученным в предыдущем пункте результатом.

Возьмём первую матрицу. разделим каждую её строку на куски размера [math]k[/math]. Для каждого куска определим номер двоичного вектора, который соответствует числам, находящимся на этом куске. Если кусок получился неравным по длине [math]k[/math](последний кусок строки), то будем считать, что в конце в нём идут не влияющие на умножение нули. Получим матрицу [math]A'_{n \times \lceil\frac{n}{k} \rceil}[/math].

Аналогично поступим с матрицей [math]B[/math], вместо строк деля столбцы. Получим матрицу [math]B'_{\lceil\frac nk\rceil\times n}[/math].

Теперь, если вместо произведения матриц [math]A[/math] и [math]B[/math] считать произведение новых матриц [math]A'[/math] и [math]B'[/math], воспользовавшись посчитанными скалярными произведениями, то каждый элемент матрицы [math]C[/math] будет получаться уже за время, пропорциональное [math]\lceil \frac nk \rceil[/math] вместо [math]n[/math], и время произведения матриц сократится с [math]O(n^3)[/math] до [math]O(n^2 \cdot\frac nk) = O(\frac{n^3}{k}) [/math].

Оценка трудоёмкости и выбор k

Оценим трудоёмкость данного алгоритма.

• Предподсчёт скалярных произведений работает за [math]O(2^{2k}k)[/math].
• Создание матриц [math]A'[/math] и [math]B'[/math] — [math]O(N^2)[/math]
• Перемножение полученных матриц — [math]O(\frac{n^3}{k})[/math]

Итого: [math]O(2^{2k}k) + O(\frac{n^3}{k})[/math]. Приведем анализ выбора числа [math]k[/math] для получения оптимальной сложности алгоритма.

В силу возрастания функции [math]f(k) = 2^{2k}k[/math] и убывания функции [math]g(k) = \frac{n^3}{k}[/math] имеем, что сложность будет оптимальна при таком значении [math]k[/math], что [math]f(k) = g(k)[/math]. Прологарифмируем обе части этого равенства:

[math]k \ln 4 + \ln k= 3 \ln n - \ln k[/math]

[math]k = \frac{3 \ln n - 2 \ln k}{\ln 4} [/math]

[math] k = 3 \log_4 n - 2 \log_4 k [/math]

В силу того, что [math] \log_4 k [/math] пренебрежительно мал по сравнению с [math] k [/math] имеем, что [math] k [/math] с точностью до константы равен [math] \log n [/math]

Таким образом, при подстановке [math]k = \log n[/math], получаем итоговую трудоёмкость [math]O(n^2 \log n) + O(\frac{n^3}{\log n}) = O(\frac{n^3}{\log n})[/math]

Код алгоритма

int n, cur;
vector <vector <int> > a, b, preculc, anew, bnew, ans;
int main() {
  freopen("input.txt", "r", stdin);
  freopen("output.txt", "w", stdout);
  cin >> n;
  a.resize(n);
  b.resize(n);
  ans.resize(n);
  // Чтение матриц
  for (int i = 0; i < n; i++) 
     for (int j = 0; j < n; j++) {
        cin >> cur;
        a[i].push_back(cur);
     }   
  for (int i = 0; i < n; i++) 
     for (int j = 0; j < n; j++) {
        cin >> cur;
        b[i].push_back(cur);
     }   
  // Предподсчёт скалярных произведений
  int k = ceil(log( (double) n));
  preculc.resize(1 << k);
  for (int i = 0; i < (1 << k); i++)
     for (int j = 0; j < (1 << k); j++) {
        int scalmul = 0;
        for (int pos = 0; pos < k; pos++)
           if (((1 << pos) & i) != 0 && ((1 << pos) & j) != 0) {
              scalmul = (scalmul + 1) % 2;
           }
        preculc[i].push_back(scalmul);
     }
  // Создание сжатых матриц
  int m = ceil(((double) n) / k);
  anew.resize(n);
  bnew.resize(m);
  for (int i = 0; i < n; i++) {
     
     for (int start = 0; start < n; start += k) {
        int cursuma = 0, cursumb = 0, curpos = start, deg = (1 << (k - 1));
        while (curpos < start + k && curpos < n) {
           cursuma += a[i][curpos] * deg;
           cursumb += b[curpos][i] * deg;
           deg /= 2;
           curpos++;
        }
        anew[i].push_back(cursuma);
        bnew[start / k].push_back(cursumb);
     }
  }
  //Перемножение полученных матриц
  for (int i = 0; i < n; i++)
     for (int j = 0; j < n; j++) {
        int curans = 0;
        for (int pos = 0; pos < m; pos++) {
           curans = (curans + preculc[anew[i][pos]][bnew[pos][j]]) % 2;
        }
        ans[i].push_back(curans);
     }
  // Вывод ответа
  for (int i = 0; i < n; i++) {
     for (int j = 0; j < n; j++) {
        cout << ans[i][j] << " ";
     }
     cout << endl;
  }
  return 0;
}