Atcoder agc039F

场上肛C肛到自闭，场后一看EF都简单过C。。。
考虑我们先给每行每列钦定一个最小值，这样每个位置能填的就是对应不小于行和列最小值的$\max$的数。然后这个显然会出问题，因为有可能某行或某列实际的最小值比钦定的大。那么考虑容斥，对于某一行或某一列，我们可以令它被容斥掉，也就是对答案带一个$-1$，并且要求填的数大于钦定的最小值。
这样就可以DP了，我们按权值从小到大考虑，令$F[i][j][k]$表示考虑了钦定最小值不超过$i$的行和列，用了$j$行和$k$列的方案数，转移枚举钦定最小值为$i+1$的没被容斥的行和列的数目及被容斥的行和列的数目，乘上容斥系数和组合数以及幂之类的一坨东西。
这样做复杂度非常高，不过容易发现没必要同时枚举那么多东西转移，我们可以依次枚举没被容斥的行数目，没被容斥的列数目，被容斥的行数目，被容斥的列数目，时间复杂度就优化到了$\mathcal{O}(nmk(n+m))$。
然后发现T爆了，因为每次转移的取模次数过于大，我的解决方案是发现转移只跟三个变量有关，于是做一个三次方的预处理，就跑快了很多。

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;

int MOD;

ll pow_mod(ll x,int k) {
  ll ans=1;
  while (k) {
  	if (k&1) ans=ans*x%MOD;
  	x=x*x%MOD;
  	k>>=1;
  }
  return ans;
}

inline void update(int &x,ll y) {
  x=(x+y)%MOD;
}

ll powd[105][10005],facd[105],facv[105];

void pre(int n,int m,int k) {
  for(int i=0;i<=k;i++) {
  	powd[i][0]=1;
  	for(int j=1;j<=n*m;j++) powd[i][j]=powd[i][j-1]*i%MOD;
  }
  facd[0]=1;
  for(int i=1;i<=max(n,m);i++) facd[i]=facd[i-1]*i%MOD;
  facv[max(n,m)]=pow_mod(facd[max(n,m)],MOD-2);
  for(int i=max(n,m)-1;i>=0;i--) facv[i]=facv[i+1]*(i+1)%MOD;
}

int f[2][105][105];
ll trans[105][105];

int main() {
  int n,m,K;
  scanf("%d%d%d%d",&n,&m,&K,&MOD);
  pre(n,m,K);
  int cur=0;
  f[0][0][0]=1;
  for(int i=1;i<=K;i++) {
  	cur^=1;
  	memset(f[cur],0,sizeof(f[cur]));
  	for(int j=0;j<=n;j++)
  	  for(int k=0;k<=m;k++) trans[j][k]=facv[j]*powd[i][j*(m-k)]%MOD*powd[K-i+1][j*k]%MOD;
  	for(int j=0;j<=n;j++)
  	for(int k=0;k<=m;k++)
  	for(int l=0;l<=j;l++)
  	    update(f[cur][j][k],f[cur^1][j-l][k]*trans[l][k]);
  	cur^=1;
  	memset(f[cur],0,sizeof(f[cur]));
  	for(int j=0;j<=m;j++)
  	  for(int k=0;k<=n;k++) trans[j][k]=facv[j]*powd[i][j*(n-k)]%MOD*powd[K-i+1][j*k]%MOD;
  	for(int j=0;j<=n;j++)
  	for(int k=0;k<=m;k++)
  	for(int l=0;l<=k;l++)
  	    update(f[cur][j][k],f[cur^1][j][k-l]*trans[l][j]);
  	cur^=1;
  	memset(f[cur],0,sizeof(f[cur]));
  	for(int j=0;j<=n;j++)
  	  for(int k=0;k<=m;k++) trans[j][k]=facv[j]*powd[i][j*(m-k)]%MOD*powd[K-i][j*k]%MOD*((j&1)?MOD-1:1)%MOD;
  	for(int j=0;j<=n;j++)
  	for(int k=0;k<=m;k++)
  	for(int l=0;l<=j;l++)
  	    update(f[cur][j][k],f[cur^1][j-l][k]*trans[l][k]);
  	cur^=1;
  	memset(f[cur],0,sizeof(f[cur]));
  	for(int j=0;j<=m;j++)
  	  for(int k=0;k<=n;k++) trans[j][k]=facv[j]*powd[i][j*(n-k)]%MOD*powd[K-i][j*k]%MOD*((j&1)?MOD-1:1)%MOD;
  	for(int j=0;j<=n;j++)
  	for(int k=0;k<=m;k++)
  	for(int l=0;l<=k;l++)
  	    update(f[cur][j][k],f[cur^1][j][k-l]*trans[l][j]);
  }
  printf("%lld\n",f[cur][n][m]*facd[n]%MOD*facd[m]%MOD);
  return 0;
}