前言

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

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

概述

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

二维分块

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

一维分块

为什么上面的方法不优呢？因为 $S⩽\omega$，所以 $\mathtt{bitset}$ 的效用发挥不出来了。$\mathtt{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}(\frac{n}{S}\cdot 2^S\cdot \frac{n}{\omega })$ 的，而 $A$ 中只需要 $\mathcal{O}(\frac{n^2}{S})$ 次行向量操作，是 $\mathcal{O}(\frac{n^3}{\omega S})$ 的，取 $S=\log n$ 有最优复杂度 $\mathcal{O}(\frac{n^3}{\omega \log n})$，且代码实现极其简单。

代码实现

$n$ 是矩阵大小，$step$ 是块的大小，已在主函数内求出。

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;
	}
};