【杂谈】01矩阵乘法

OneInDark

已于 2022-08-05 19:31:00 修改

阅读量962

点赞数 4

分类专栏： C++ # 状态压缩分块文章标签：太没用了

于 2021-11-09 14:33:25 首次发布

本文链接：https://blog.csdn.net/qq_42101694/article/details/121227383

版权

C++ 同时被 3 个专栏收录

547 篇文章 3 订阅

订阅专栏

状态压缩

23 篇文章 0 订阅

订阅专栏

分块

13 篇文章 0 订阅

订阅专栏

本文探讨了优化01矩阵乘法的一种一维分块方法，通过在模2意义下进行计算，利用bitset优化，实现了比传统二维分块更优的复杂度O(n³/ωlogn)。代码实现简洁，尤其适用于大块大小。通过在kk方向进行分块，减少了计算量，达到优于二维分块的效果。

摘要由CSDN通过智能技术生成

前言

前几天看了一下 $2017$ 年论文《非常规大小分块算法初探》，里面有一节讲到了 $01$ 矩阵乘法，介绍了一种 $\mathcal O({n^3\over\omega\sqrt{\log n}})$ 的分块做法。然而其代码实现难度高，复杂度优化效果小，我曾向 $\sf Rainybunny$ 说，不会真的有人写这玩意儿骗分吧！

结果我还是写了分块 $01$ 矩阵乘法来骗分。但是并不是上面的方法。是我无意中想到的一个分块方法。为什么要写出来呢？可能也是为了留作纪念吧……

概述

这里的 $01$ 矩阵是，在模 $2$ 意义下（在域 $\Bbb{GF}(2)$ 内）进行所有操作，包括答案矩阵也要对 $2$ 取模。

二维分块

这是论文中的方法，先简单讲讲。不妨设要计算 $A\times B$ ，将两个矩形都划分成若干 $S\times S$ 大小的子矩形。花费 $\mathcal O(4^{S^2}{S^3\over \omega})$ 的时间，预处理每一种不同的情况对应的结果。然后用 $\mathcal O[({n\over S})^3{S^2\over\omega}]=\mathcal O({n^3\over S\omega})$ 的复杂度合并结果。取 $S=\sqrt{{\log_4 n}\over 2}$ 时，有最优复杂度 $\mathcal O({n^3\over{\omega \sqrt{\log n}}})$ 。

一维分块

为什么上面的方法不优呢？因为 $S\leqslant\omega$ ，所以 $\tt bitset$ 的效用发挥不出来了。 $\tt bitset$ 才是最强的分块者！

考虑 $\omega$ 和 $S$ 的分块效果同时存在。回到原式 $C_{i,j}=\sum_k A_{i,k}B_{k,j}$ 上，不妨设 $i$ 是未分块的一维。能不能对 $j$ 做 $S$ 分块呢？相当于矩阵 $A$ 的某个行向量，要点乘 $S$ 个列向量。似乎并没有冗余计算。估计行不通。

所以考虑对 $k$ 做 $S$ 分块。也就是 $A$ 的行向量的一部分，长度仅为 $S$ ，代表 $B$ 的某些行向量是否被选择，然后异或到第 $i$ 行上。长度仅为 $S$ ，这提示我们 $\text{Method of four Russians}$ ，即讨论 $2^S$ 种不同的 $A_{i,k}\;(k\in[l,r])$ 的情况，并直接计算出对第 $i$ 行的贡献。

那么 $B$ 中每 $S$ 行就要暴力计算一次，这是 $\mathcal O({n\over S}\cdot 2^S\cdot{n\over\omega})$ 的，而 $A$ 中只需要 $\mathcal O({n^2\over S})$ 次行向量操作，是 $\mathcal O({n^3\over \omega S})$ 的，取 $S=\log n$ 有最优复杂度 $\mathcal O({n^3\over \omega\log n})$ ，且代码实现极其简单。

代码实现

$n$ 是矩阵大小， $s t e p$ 是块的大小，已在主函数内求出。

const int MAXN = 2000;
const int LOGN = 10;
bitset<MAXN> dp[1<<LOGN];
int lowbit[1<<LOGN], n, step;
struct Matrix{
	bitset<MAXN> a[MAXN];
	inline void clear(){
		for(int i=0; i!=n; ++i)
			a[i].reset();
	}
	Matrix operator * (const Matrix &b) const {
		Matrix c; c.clear();
		for(int l=0,r=step; l!=n; l=r,r+=step){
			if(__builtin_expect(r > n, false))
				r = n; // not to exceed n
			for(int S=1; S!=(1<<step); ++S){
				const int &id = lowbit[S];
				dp[S] = dp[S^(1<<id)]^b.a[l+id];
			}
			for(int i=0,now=0; i!=n; ++i,now=0){
				for(int j=l; j!=r; ++j)
					now ^= (a[i].test(j))<<(j-l);
				c.a[i] ^= dp[now]; // find in table
			}
		}
		return c;
	}
};