不基于 FKT Algorithm 的网格图完美匹配计数做法

参考:http://www.algonotes.com/en/counting-matchings/ 和四川省选 D 官方题解。

下午研究 FKT Algorithm 的时候发现了这篇博客,弄明白后觉得非常巧妙,而且实现极为容易!

仍然有一些只会感性理解不会严格证明的结论…

题目描述

给定一个 $n \times m$ 的棋盘,求有多少种在棋盘上放置 $1 \times 2$ 的多米诺骨牌使棋盘密铺的方法。同时,对于每个位置 $(i,j)$ ,分别给出在 $(i,j)\to (i,j+1)$ 和 $(i,j) \to (i + 1,j)$ 的放置方案数 $R_{i,j},D_{i,j}$ ,保证不存在 $i,j$ 使得 $R_{i,j}D_{i,j} = 0$ 。求出答案对 $998244353$ 取模的值。

保证 $2 \le n,m \le 300$ 。

二分图完美匹配计数

原问题等价于一个 $n \times m$ 的网格图上求完美匹配的个数。

可以通过将 $(1,1)$ 设为白点,与之相邻的点设为黑点再不断进行染色的方法将一个棋盘黑白染色。然后放置的骨牌一定同时覆盖了一个黑点和一个白点,也就是说网格图上可以选择的边的两端颜色不同,故可以转化成二分图求完美匹配方案数的问题。很显然,如果 $nm$ 为奇数则答案一定为 $0$ ,否则二分图左右分别有 $\frac{nm} 2$ 个点。

二分图求完美匹配方案数是一个经典的 Hard Problem 。我们考虑一个矩阵 $A$ 满足如果左边的 $i$ 点到右边的 $j$ 点有边则置 $A_{i,j} = w_{i,j}$ ,这里的 $w_{i,j}$ 就是这个位置放置骨牌的方案数,那么这个矩阵的积和式就是完美匹配的方案数(下面有两个例子)。而求积和式是没有多项式算法的(曾经写过 积和式模 4),所以这个做法难以解决问题。

当然,这样可以得到答案模 $2$ 的值,因为此时行列式和积和式等价。

网格图完美匹配计数

我们的目标是使用行列式而不是积和式来解决问题,因为前者是可以在多项式时间计算的。FKT Algorithm 会把两个任意完美匹配拼起来得到一个偶环剖分从而通过定向使得答案为 $\sqrt{\det G}$ ,但是网格图中有更加巧妙地做法。

首先,根据行列式的定义,对于一个奇排列,其贡献为 $-\prod A_{i,\pi_i}$ ,而偶排列则为 $\prod A_{i,\pi _i}$ 。同时,奇排列中会有奇数个偶环,而偶排列中有偶数个偶环,故我们希望对于一个有奇数个偶环的排列,通过调整边权,使得其另外获得 $-1$ 的贡献,而有偶数个偶环的排列不发生变化。也就是说,如果能够让每个偶环都另外造成 $-1$ 的贡献,那么就可以用行列式求出答案了。

接下来考虑怎么在网格图上表示“排列中的偶环”。现在我们考虑一个偶环 $\pi = [3,1,4,2]$ :

也就是说,如果我们像红边一样将左右对应的点全部连起来,那么一个排列中的偶环对应的匹配就会扩展成原图上的偶环,当然,这里也可能出现长为 $2$ 的偶环。

接下来,对于一个选定的完美匹配,在网格图上我们将匹配中的边和二分图中对应点之间的红边均加粗,然后会得到一个偶环剖分。这里如果 $m$ 不为偶数可以先转秩一下,然后让水平相邻点两两配对成为二分图中的对应点。

其中 $(2,2),(5,5)$ 是长度为 $2$ 的偶环。

然后,我们将所有竖直边上的权值从 $w_{u,v}$ 变成 $w_{u,v}i$ ,这里 $i$ 为虚数单位,直接求行列式后得到的就是答案。

我们将一个原图上的选中的偶环里的边分为三类:

  • $e_H$ 水平边,同时存在于完美匹配中。
  • $e_B$ 水平边,但是不存在于完美匹配中,也就是加入的连接左右对应点的红边,也是上图中 $(1,1),(4,4)$ 这种边。
  • $e_V$ 竖直边,它一定在完美匹配中。

那么我们要说明的就是 $i^ = -1$ ,因为在原来的权值上新加入的权值就是竖直边上 $i$ 的乘积,而由于网格图上的环一定是个偶环,所以新增加的权值肯定是 $-1$ 。

首先,$e_V,e_H+e_B$ 一定是偶数。考虑按任意顺序从一个点出发走完整个环,向下一次就一定会向上一次,左右同理。

然后,对于两个行标号奇偶性不同的 $e_B$ 边,在遍历过程中经过这两条边的方向一定不同。不妨假设两个 $e_B$ 边按从靠左点到靠右点分别是 $a \to b,c \to d$ ,那么在从 $a$ 到达 $c$ 的过程中列的标号的奇偶性是不会改变的,而行的奇偶性改变了,那么所经过的边数为奇数,而偶环中一定是按照 $e_B,e_H/e_V$ 交替来经历所有边的,所以一定不可能以 $c \to d$ 的顺序经过这条 $e_B$ 边。

又显然有经过最上面的行和最下面的行的时候方向是相反的,所以上下行数差为奇数,所以 $\frac{e_V}{2}$ 为奇数。

这部分的证明写的不是很严谨,因为确实不太会严格证明这个

因此我们在这个图上求行列式即可。$i$ 可以真的看成复数做计算,但是事实上模 $998244353$ 时 $-1$ 有二次剩余,因此可以直接用 $86583718$ 做为 $i$ 。

这样就有了一个 $O(n^6)$ 的做法。

优化高斯消元

这部分利用了题解中所说的主元法,它能应用依赖于所有边权都不为 $0$ 。

其实这个矩阵是很稀疏的, $d = (x-1)m’+y$ 相关的点只有 $d,d-m’,d+m’$ 以及 $d-1,d+1$ 中的一个(取决于当前行为奇还是偶)。我们考虑把每一行都消成只有 $[1,m’]$ 以及 $d+m’$ 这些位置有值。具体来说,对于 $t=d,d-m’,d-1/d+1$ ,我们一定已经消元了第 $t-m’$ 行,在这一行中会剩下 $t$ 这个位置的值(由于边权非 $0$ ,这个值一定不是 $0$ ),用这个值来消掉对应的 $t$ 即可。而对于前 $m’$ 行,直接存下 $d,d-m’,d-1/d+1$ 即可。

这个过程中,计算 $d$ 行时最多只会用到 $[d-2m,d]$ 的行,因此只需要存 $O(m)$ 个行即可,空间只需要 $O(m^2)$ 。

然后对于 $[m+1,nm]$ 这些列,每列一定只有一个值,所以行列式的时候只能选对应值,直接乘进答案即可,然后就可以删除掉对应的行,最后只有一个大小为 $m’$ 的矩阵参与高斯消元。

需要注意的是,之前那些被固定的列中的值仍然会对排列的奇偶造成影响,所以需要乘上 $(-1)^{(n-1)m’^2}$ 。

时间复杂度 $O(m^3 + nm)$ ,空间复杂度 $O(m^2 + nm)$ 。由于出题人不卡空间,实现直接开了 $O(m^3)$ 的空间。

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
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
typedef long long ll;
const int P = 998244353;
const int I = 86583718;
const int MAXN = 306;
int n , m , mm;

int qm( int x ) {
return x < P ? x : x - P;
}

int Pow( int x , int a = P - 2 ) {
int ret = 1;
while( a ) {
if( a & 1 ) ret = ret * 1ll * x % P;
x = x * 1ll * x % P , a >>= 1;
}
return ret;
}

int tA[MAXN][MAXN] , tB[MAXN][MAXN] , A[MAXN][MAXN] , B[MAXN][MAXN];

struct vec {
int A[MAXN];
vec( ) { mem( A ); }
int& operator [] ( const int x ) { return A[x]; }
vec operator + ( vec t ) {
vec ret;
rep( i , 1 , mm ) ret[i] = qm( A[i] + t[i] );
return ret;
}
vec operator - ( vec t ) {
vec ret;
rep( i , 1 , mm ) ret[i] = qm( A[i] + P - t[i] );
return ret;
}
vec operator * ( int x ) {
vec ret;
rep( i , 1 , mm ) ret[i] = A[i] * 1ll * x % P;
return ret;
}
} V[MAXN * MAXN] ;

int gid( int x , int y ) {
return ( x - 1 ) * mm + ( y );
}

int vid[MAXN * MAXN];
int G[MAXN][MAXN];

int det( int n ) {
int ans = 1;
rep( i , 1 , n ) rep( j , 1 , n ) G[i][j] = ( G[i][j] + P ) % P;
rep( i , 1 , n ) {
if( !G[i][i] ) {
rep( j , i + 1 , n ) {
if( G[j][i] ) {
swap( G[j] , G[i] ) , ans = P - ans; break;
}
}
}
ans = ans * 1ll * G[i][i] % P;
rep( j , i + 1 , n ) if( G[j][i] ) {
int kk = G[j][i] * 1ll * Pow( G[i][i] , P - 2 ) % P;
rep( k , i , n )
G[j][k] = ( G[j][k] + P - 1ll * kk * G[i][k] % P ) % P;
}
}
return ans;
}

void solve() {
cin >> n >> m;
rep( i , 1 , n ) rep( j , 1 , m - 1 ) scanf("%d",tA[i] + j);
rep( i , 1 , n - 1 ) rep( j , 1 , m ) scanf("%d",tB[i] + j);
if( ( n & 1 ) && ( m & 1 ) ) {
puts("0"); return;
}
if( m & 1 ) {
swap( n , m );
rep( i , 1 , n ) rep( j , 1 , m - 1 ) A[i][j] = tB[j][i];
rep( i , 1 , n - 1 ) rep( j , 1 , m ) B[i][j] = tA[j][i];
} else {
memcpy( A , tA , sizeof A ) , memcpy( B , tB , sizeof B );
}
mm = m / 2;
rep( i , 1 , mm ) {
vid[i] = B[1][i * 2 - 1] * 1ll * I % P;
V[i][i] = A[1][i * 2 - 1];
if( i != 1 ) V[i][i - 1] = A[1][i * 2 - 2];
}
rep( i , 2 , n ) rep( j , 1 , mm ) {
int d = gid( i , j ) , tj = ( j * 2 - ( i & 1 ) );
V[d] = V[d] - V[d - mm] * ( A[i][( i & 1 ) ? tj : tj - 1] * 1ll * Pow( vid[d - mm] ) % P );
if( i & 1 ) {
if( tj != 1 )
V[d] = V[d] - V[d - 1 - mm] * ( A[i][tj - 1] * 1ll * Pow( vid[d - 1 - mm] ) % P );
} else {
if( tj != m )
V[d] = V[d] - V[d + 1 - mm] * ( A[i][tj] * 1ll * Pow( vid[d + 1 - mm] ) % P );
}
if( i > 2 )
V[d] = V[d] - V[d - mm - mm] * ( B[i - 1][tj] * 1ll * I % P * Pow( vid[d - mm - mm] ) % P );
else
V[d][d - mm] = ( V[d][d - mm] + B[i - 1][tj] * 1ll * I ) % P;
if( i != n ) vid[d] = B[i][tj] * 1ll * I % P;
}
// rep( i , 1 , n * mm ) {
// rep( j , 1 , mm ) printf("%d ",V[i][j]); puts("");
// printf("%d\n",vid[i]); puts("");
// }
int res = 1;
rep( i , 1 , ( n - 1 ) * mm ) res = res * 1ll * vid[i] % P;
rep( i , ( n - 1 ) * mm + 1 , n * mm ) rep( j , 1 , mm ) G[i - ( n - 1 ) * mm][j] = V[i][j];
res = res * 1ll * det( mm ) % P;
if( ( n * mm - mm ) * mm & 1 ) printf("%d\n",P - res);
else printf("%d\n",res);
}

signed main() {
// freopen("in","r",stdin);
// freopen("out","w",stdout);
// int T;cin >> T;while( T-- ) solve();
solve();
}


一个可能更棒的解释:

如果 $n,m$ 均为奇数,答案为 $0$ 。

如果 $m \equiv 1 \pmod 2$ 则转秩一下,使得 $m$ 为偶数。然后对网格图进行黑白染色(二分图染色),匹配的边的两端颜色一定不同。

建立矩阵 $G$ ,满足 $G_{i,j}$ 为 $1$ 当且仅当左部的 $i$ 到右部的 $j$ 有变。二分图完美匹配的数量即求有多少排列满足矩阵的 $(i,p_i)$ 位置为 $1$ ,这就是矩阵的积和式,难以计算。考虑怎么把积和式转化成行列式,后者可以在多项式复杂度得出值。

事实上,只需要对于所有在网格图中竖着的边的权值由 $1$ 修改为 $i = \sqrt{-1}$ 即可直接求行列式来得到积和式的值。具体证明如下。

首先,我们钦定所有左右部对应的点之间的边被“选中”,如图:

image.png

然后在网格图中随意加入一个完美匹配,会发现完美匹配与被“选中”的边在一起会形成若干个环,并且这里的环其实就是排列中通过 $i \to p_i$ 形成的环再加入“选中”的边得到的。同时,得到的偶环一定满足环上的边由来自完美匹配的边和来自“选中”的边交替形成。

那么,如果我们定义一个偶环中:

  • $e_B$ 表示环里被选中的边的数量;
  • $e_H$ 表示环里完美匹配中的横边数量;
  • $e_V$ 表示环里竖直边的数量。

那么由于交替出现,有 $e_V + e_H = e_B$ ,同时由于从一个点出发最终会回到自己,有 $e_V,e_B+e_H$ 均为偶数,因为两者分别能对列数、行数的奇偶性发生改变。我们考虑对于相邻边,如果是 $e_B \to e_H$ 那么会对当前列数改变 $2$ ,但是 $e_B \to e_V$ 永远只能改变当前列数的奇偶性,类似于 $\text{xor}\ 1$ ,因此 $e_B \to e_H$ 这种操作的个数也是偶数,只有这样才能让列数除以 $2$ 的值回到原点。所以 $e_H$ 也是偶数,故 $e_B,e_H,e_V$ 均为偶数。所以说,对于求积和式时排列中的任何一个环,环的大小一定是偶数。

回忆一下行列式的定义,我们需要让每个排列中的偶环(即由 $i \to p_i$ 连接而形成的偶环)通过矩阵中的值的乘积另外造成 $-1$ 的贡献。而我们证明了每个有值排列的环都是偶环,也就是说只需要让每个排列中的环都造成 $-1$ 的贡献即可。

现在我们证明偶环中 $\frac{e_V}{2}$ 一定为奇数。

首先,如果两个行的标号的奇偶性不同,那么在经过这两个行中的“选中”的方向一定相反。因为交替经过偶环中的边的时候一定是从二分图左边的 $i$ 到右边的 $p_i$ 再经过 $R_{p_i} \to L_{p_i}$ (或者反过来)回到左边,因此所有被选中边都是从左部到右部经过,而由图可知如果两行行号奇偶性不同,则“选中”边的左右关系相反。

然后,考虑从任意一个行号最小的点中列号最小的点 $u$ 出发,绕着偶环走一圈回到该点的过程。我们把过程分为从 $u$ 到行号最小的列号最大的点、环上行号最大的某点,再从该点回到行号最小点的过程。由于得到的是一个环,经过最上方和最下方时的方向肯定是不同的,因此两行号差为奇数。而一个垂直边会对行号造成 $1$ 的贡献,因此两过程中的 $e_V$ 均为奇数,因此 $i^{\frac{e_V}{2}} = -1$ ,故每个排列中的偶环,也即原图中的偶环,会对行列式的值造成 $-1$ 的贡献,从而抵消掉行列式计算时的 $-1$ 贡献,得到积和式的值。

于是就可以把矩阵求出后直接 $O(n^3m^3)$ 得到行列式的值。

考虑利用主元法去求行列式的值。我们尝试把每一行都消到只有 $(1,1),\ldots,(1,m’)$ 和 $(i+1,j)$ 这 $m’ + 1$ 个变量。

考虑到对于 $(i,j)$ 位置的点,与其有边的只有 $(i-1,j)(i,j-1),(i,j+1),(i+1,j)$ 的位置,也就是说这行只有这四列有值。于是我们考虑把这一行通过之前的行消元,使得其只有 $(1,1),\ldots,(1,m’),(i+1,j)$ 这些列有值,然后用这行去消元 $(i+1,j-1),(i+1,j+1),(i+2,j)$ 对应的行。

对于 $j > m’$ 的所有列 $j$ ,其中一定只有一个行对应有值,因此在求行列式的时候这些位置必选。剩下的就只有 $(n-m’,1)$ 到 $(n,m’)$ 这个矩阵有值了,直接进行 $O(m^3)$ 的高斯消元即可。记得最后乘上之前选的位置对逆序对个数的贡献。

复杂度 $O(m^3)$ 。

注意到 $-1$ 在模 $998244353$ 存在二次剩余,直接取 $i = 86583718$ 即可。

文章作者: yijan
文章链接: https://yijan.co/domino/
版权声明: 本博客所有文章除特别声明外,均采用 CC BY-NC-SA 4.0 许可协议。转载请注明来自 Yijan's Blog