TC13986 SubRectangles加强版

把 \(A\) 的范围由 \(4\) 加强到了 \(24\)。

Problem

给定 \(H\) 行 \(W\) 列的矩阵 \(a\)，每个位置可以填 \(0\) 或 \(1\)，一个矩阵合法当且仅当其任意一个 \(A\times B\) 的子矩阵的和相同，求合法矩阵数。

\(H,W\le 10^9,A\le 24,B\le 4\)

Solution

对于一个点 \((i,j)\)，如果 \((i-A,j),(i,j-B),(i-A,j-B)\) 已经确定，则这个格子唯一确定。证明：

上图为 \(4\) 个 \(A\times B\) 的矩形，于是有 \(A+E+H=E+B+F=H+C+G=F+G+D\)，从而有 \(A+E=C+G\)，\(D+G=B+E\)，化简后有 \(D=B+C-A\)。

那么只需要填完前 \(A\) 行和前 \(B\) 列就可以唯一确定整个矩形。

但是会有不合法的情况，也就是 \(a_{i-A,j}=a_{i,j-B}=1-a_{i-A,j-B}\) 时， \(a_{i,j}\) 为 \(-1\) 或 \(2\)。

为了进一步分析性质，将 \((i,j)\) 按 \(i\bmod A\) 和 \(j\bmod B\) 分组后拆成 \(A\times B\) 个矩阵，每个重标号为矩阵 \(b\)。那么上述限制变为 \(b_{i,j}=b_{i-1,j}+b_{i,j-1}-b_{i-1,j-1}\)，\(b\) 矩阵合法的条件是 \(b_{i-1,j}= b_{i-1,j-1}\) 或 \(b_{i,j-1}=b_{i-1,j-1}\)。

引理：若 \(i>0,j>0\)，则 \(b_{i,j}=b_{i,0}+b_{0,j}-b_{0,0}\)。

证明：按 \(i,j\) 从小到大归纳，对于 \((i',j')\ne (i,j)\)，若 \(i'\le i\) 且 \(j'\le j\) 则结论成立，要证 \((i,j)\) 成立。

1. 若 \(i>1,j>1\)，则 \(b_{i,j}=b_{i-1,j}+b_{i,j-1}-b_{i-1,j-1}=b_{i,0}+b_{0,j}-b_{0,0}\)
2. 若 \(i=1,j=1\)，显然成立。
3. 若 \(i=1\)，则 \(b_{i,j}=b_{0,j}+b_{i,j-1}-b_{0,j-1}=b_{i,0}+b_{0,j}-b_{0,0}\)
4. 若 \(j=1\)，与 \(i=1\) 同理。

于是矩阵合法转化为任意 \(i,j>0\)，有 \(b_{i,0}= b_{0,0}\) 或 \(b_{0,j}=b_{0,0}\)，从而推出要不然第 \(0\) 行全等于 \(b_{0,0}\)，要不然第 \(0\) 列全等于 \(b_{0,0}\)。

回到原问题，考虑前 \(A\times B\) 个格子，每个格子有三种选择，分别是：行相同，列相同，都相同，记作 \(1,2,3\)。

考虑如果确定了每个格子的选择，如何计算方案数。

因为有 \(3\) 类格子，所以需要容斥，容斥系数为 \((-1)^{x}\)，其中 \(x\) 是 \(3\) 类格子的个数。

现在的限制是，第 \(i\) 行前 \(B\) 个的和要与第 \(i+A\) 行前 \(B\) 个的和相同，列也有同样的限制。

因为 \(3\) 类格子不影响限制，所以方案数先乘 \(2^x\)，然后就可以不管 \(3\) 类格子。

考虑行列相互独立，拆开来算，以行的限制计算为例：

记与第 \(i\) 行\(\bmod A\) 相同的行有 \(w\) 个，这个可以 \(O(1)\) 算。记第 \(i\) 行有 \(j\) 个 \(1\) 类格子，\(2\) 类格子和 \(3\) 类格子由于所有行在这列一样，所以不用管。

枚举 \(j\) 个 \(1\) 类格子有 \(k\) 个选 \(1\)，则每行独立且都要选 \(k\) 个，方案数 \(\dbinom{j}{k}^{w}\)，所以总贡献为 \(\sum\limits_{k=0}^{j}\dbinom{j}{k}^{w}\)

每个 \(\bmod A\) 不同的行之间独立，所以贡献乘起来就行。

对于列的计算是几乎一致的，于是我们有 \(O(3^{AB}Poly)\) 的做法。

下面这份代码有很多东西没预处理，所以无法通过 \(A,B\le 4\)，但这不重要。

#include<bits/stdc++.h>
using namespace std;
#define int long long
int H,W,A,B,ans,a[30][5],C[50][50];
const int mod=998244353;
int q_pow(int x,int y){
	if(!y) return 1;
	int z=q_pow(x,y>>1);
	if(y&1) return z*z%mod*x%mod;
	else return z*z%mod;
}
void solve(){
	int cnt=0,sum=1;
	for(int i=1;i<=A;i++){
		for(int j=1;j<=B;j++){
			if(a[i][j]==3){
				cnt++;
				sum=sum*2%mod;
			}
		}
	}
	for(int i=1;i<=A;i++){
		int x=0,w=(H-i)/A+1,y=0;
		for(int j=1;j<=B;j++) if(a[i][j]==1) x++;
		for(int j=0;j<=x;j++){
			y=(y+q_pow(C[x][j],w))%mod;
		}
		sum=sum*y%mod;
	}
	for(int i=1;i<=B;i++){
		int x=0,w=(W-i)/B+1,y=0;
		for(int j=1;j<=A;j++) if(a[j][i]==2) x++;
		for(int j=0;j<=x;j++){
			y=(y+q_pow(C[x][j],w))%mod;
		}
		sum=sum*y%mod;
	}
	ans=(ans+sum*(cnt&1?mod-1:1))%mod;
}
void dfs(int x,int y){
	if(x>A){
		solve();
		return;
	}
	if(y>B){
		dfs(x+1,1);
		return;
	}
	a[x][y]=1;
	dfs(x,y+1);
	a[x][y]=2;
	dfs(x,y+1);
	a[x][y]=3;
	dfs(x,y+1); 
}
signed main(){
	for(int i=0;i<=30;i++){
		C[i][0]=1;
		for(int j=1;j<=i;j++) C[i][j]=(C[i-1][j-1]+C[i-1][j])%mod;
	}
	cin>>H>>W>>A>>B;
	if(A>H || B>W){
		cout<<q_pow(2,H*W);
		return 0;
	}
	dfs(1,1);
	cout<<ans;
	return 0;
}

考虑如何优化，这一步比较自然，观察上面的贡献计算，可以发现我们只关心每行有几个 \(1\) 类格子，每列有几个 \(2\) 类格子，以及每个 \(3\) 类格子带来 \(-2\) 的贡献。

因为 \(B\) 很小，考虑从上往下 dp，记 \(f_{i,S}\) 表示当前在第 \(i\) 行，每列 \(2\) 类格子数量的状态为 \(S\) 的方案数。

枚举下一行哪些列选 \(2\) 类格子，下一行选了几个 \(1\) 类格子后即可转移，由于转移种类很少，预处理后容易做到 \(O((A+1)^{B}A2^B)\) 的复杂度。

这东西也容易用轮廓线 dp 做到 \(O((A+1)^{B}AB^2)\) 的复杂度，只需要在状态中多记一个这一行选了几个 \(1\) 类格子即可。

下面是 \(O((A+1)^{B}A2^B)\) 的代码。

#include<bits/stdc++.h>
using namespace std;
#define int long long
int H,W,A,B,ans,C[50][50],val[50][50],val2[50][50],f[30][400010];
const int mod=998244353;
int q_pow(int x,int y){
	if(!y) return 1;
	int z=q_pow(x,y>>1);
	if(y&1) return z*z%mod*x%mod;
	else return z*z%mod;
}
signed main(){
	for(int i=0;i<=30;i++){
		C[i][0]=1;
		for(int j=1;j<=i;j++) C[i][j]=(C[i-1][j-1]+C[i-1][j])%mod;
	}
	cin>>H>>W>>A>>B;
	if(A>H || B>W){
		cout<<q_pow(2,H*W);
		return 0;
	}
	int V=25*25*25*25;
	f[1][0]=1;
	for(int i=1;i<=A;i++){
		for(int j=0;j<=B;j++){
			int w=(H-i)/A+1;
			for(int x=0;x<=j;x++){
				int y=0;
				for(int k=0;k<=x;k++) y=(y+q_pow(C[x][k],w))%mod;
				y=y*q_pow(mod-2,j-x)%mod*C[j][x]%mod;
				val[i][j]=(val[i][j]+y)%mod;
			}
		}
	}
	for(int i=1;i<=B;i++){
		int w=(W-i)/B+1;
		for(int x=0;x<=A;x++) for(int j=0;j<=x;j++) val2[i][x]=(val2[i][x]+q_pow(C[x][j],w))%mod;
	}
	for(int i=1;i<=A;i++){
		for(int j=0;j<=V;j++){
			if(!f[i][j]) continue;
			for(int S=0;S<(1<<B);S++){
				int nj=j,cnt=B-__builtin_popcount(S);
				if(S&1) nj+=1;
				if(S&2) nj+=25;
				if(S&4) nj+=25*25;
				if(S&8) nj+=25*25*25;
				f[i+1][nj]=(f[i+1][nj]+f[i][j]*val[i][cnt])%mod;
			}
		}
	}
	for(int j=0;j<=V;j++){
		int k=j,sum=1;
		for(int i=1;i<=B;i++){
			sum=sum*val2[i][k%25]%mod;
			k/=25;
		}
		ans=(ans+f[A+1][j]*sum)%mod;
	}
	cout<<ans;
	return 0;
}