不基于 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();
}

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