很厉害的一道 Ad-hoc 题！假定我们填写的矩阵为 \(A\)。
直接填写 \(A\) 计算贡献是基本做不到的，考虑挖掘一些神奇的东西。
问题转化：考虑 \(\prod\limits_{i=1}^n \prod \limits_{j=1}^m f(i,j)\) 相当于构造另一个 \(B\) 矩阵，满足 \(B_{i,j}\le \min(A_{i,1...m},A_{1...n,j})\)，问矩阵对 \((A,B)\) 的个数。
继续挖掘条件，发现 \(\forall (i,j), \ B_{i,j} \le \min A_{i,1...m} \Leftrightarrow \forall (i,j) , \ A_{i,j} \ge \max B_{i,1...m}\)，列同理。
所以我们两边一起跑，设 \(x_i=\max B_{i,1...m},\ y_j=\min A_{1...n,j}\)。
而矩阵对合法的充要条件转化为 \(\forall (i,j) , \ B_{i,j} \le \min(x_i,y_j) \ \wedge \ A_{i,j} \ge \max(x_i,y_j)\)，且 \(B_{i,1...m}\) 至少一个等于 \(x_i\)，\(A_{1...n,j}\) 至少一个等于 \(y_j\)。
依次填写 \(x_{1...n},y_{1...m}\)，设 \(f[k,i,j]\) 表示填写了 \(\le k\) 的数，其中 \(x_{1...n}\) 填写了 \(x\) 个，其中 \(y_{1...m}\) 填写了 \(y\) 个。
考虑在填写 \(x_{1...n},y_{1...m}\) 时一起填写矩阵 \(A,B\)。对于 \(B_{i,j}\)，当填写了 \(x_i,y_j\) 中较小者后计算贡献；对于 \(A_{i,j}\)，当填写了 \(x_i,y_j\) 中较大者后计算贡献。
考虑从 \(f[i-1,x,y]\) 转移过来，枚举 \(p,q\) 表示 \(x_{1...n},y_{1...m}\) 中等于 \(=k\) 的数的个数。
-
对于 \(x_s=k,\ y_t<k\)，那么 \(A_{s,t}\) 贡献为 \(k-i+1\)。
-
对于 \(x_s=k,\ y_t\ge k\)，那么 \(B_{s,t}\) 贡献为 \(i\)，并且对于每个 \(s\) 都至少一个 \(B_{s,t}\) 等于 \(i\)。
-
对于 \(x_s\le k,\ y_t=k\)，那么 \(A_{s,t}\) 贡献为 \(k-i+1\)，并且对于每个 \(t\) 都至少一个 \(A_{s,t}\) 等于 \(i\)。
-
对于 \(x_s>k,\ y_t=k\)，那么 \(B_{s,t}\) 贡献为 \(i\)。
枚举的 \(p,q\) 之间显然任何关系，可以分步转移。
预处理一些必要的系数，时间复杂度 \(O(knm(n+m))\)。
启示：
-
问题转化，深入思考
-
等价条件转化，使统计变得更方便
-
利用条件，精准刻画
点击查看代码
#include<bits/stdc++.h>
#define ll long long
#define fi first
#define se second
#define mkp make_pair
#define pir pair<ll,ll>
#define ull unsigned ll
#define pb push_back
#define i128 __int128
using namespace std;
const ll maxn=110;
ll n, m, t, mod;
ll power(ll a,ll b=mod-2){
ll s=1;
while(b){
if(b&1) s=s*a%mod;
a=a*a%mod; b>>=1;
} return s;
}
ll f[maxn][maxn], g[maxn][maxn], r[maxn][maxn], c[maxn][maxn], pw[maxn][maxn];
ll P_pw[maxn][maxn][maxn], P_r[maxn][maxn][maxn];
int main(){
scanf("%lld%lld%lld%lld",&n,&m,&t,&mod);
c[0][0]=1;
for(ll i=1;i<=n||i<=m;i++){
c[i][0]=1;
for(ll j=1;j<=i;j++)
c[i][j]=(c[i-1][j]+c[i-1][j-1])%mod;
}
for(ll i=1;i<=t;i++) pw[0][i]=1;
for(ll i=1;i<=n||i<=m;i++)
for(ll j=1;j<=t;j++) pw[i][j]=pw[i-1][j]*j%mod;
for(ll i=1;i<=n||i<=m;i++)
for(ll j=1;j<=t;j++)
r[i][j]=(power(j,i)-power(j-1,i)+mod)%mod;
for(ll i=0;i<=n||i<=m;i++)
for(ll j=1;j<=t;j++){
P_pw[i][j][0]=P_r[i][j][0]=1;
for(ll k=1;k<=n||k<=m;k++)
P_pw[i][j][k]=P_pw[i][j][k-1]*pw[i][j]%mod,
P_r[i][j][k]=P_r[i][j][k-1]*r[i][j]%mod;
}
f[0][0]=1;
for(ll k=1;k<=t;k++){
memset(g,0,sizeof g);
for(ll i=0;i<=n;i++)
for(ll j=0;j<=m;j++){
for(ll x=0;i+x<=n;x++){
g[i+x][j]=(g[i+x][j]+f[i][j]*P_r[m-j][k][x]%mod*
P_pw[j][t-k+1][x]%mod*c[n-i][x])%mod;
}
}
memset(f,0,sizeof f);
for(ll i=0;i<=n;i++)
for(ll j=0;j<=m;j++){
for(ll y=0;j+y<=m;y++){
f[i][j+y]=(f[i][j+y]+g[i][j]*P_pw[n-i][k][y]%mod*
P_r[i][t-k+1][y]%mod*c[m-j][y])%mod;
}
}
} printf("%lld",f[n][m]);
return 0;
}