把 \(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)\) 成立。
- 若 \(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}\)
- 若 \(i=1,j=1\)，显然成立。
- 若 \(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}\)
- 若 \(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;
}