Link
先考虑将问题进行转化，我们要对满足如下条件的大小为(n imes m)的矩阵(A,B)进行计数：
对任意(A_{i,j}(iin[1,n],jin[1,m]))，(B)中(B_{i,j})所在行和列的每一个元素都不小于(A_{i,j})。
也就是对于任意一个矩阵(B)，它在题意中的权值等价于合法的(A)矩阵的数目。
不难发现这个限制条件与下述限制条件等价：
(forall iin[1,n],maxlimits_{j=1}^m(A_{i,j})leminlimits_{j=1}^m(B_{i,j}))
(forall jin[1,m],maxlimits_{i=1}^n(A_{i,j})leminlimits_{i=1}^n(B_{i,j}))
设(f_{h,i,j})表示(B)中有(i)列，(A)中有(i)行中的数都(le k)的方案数。
如果确定(A)中某行最大值为(h+1)，则(B)中对应行确定了列最小值的位置可以填(ge h+1)的数；没有确定列最小值的位置可以填(le h+1)的数，且至少有一个位置填(h+1)。
如果确定(B)中某列最小值为(h+1)，则(A)中对应列确定了行最大值的位置可以填(ge h+1)的数，且至少有一个位置填(h+1)；没有确定行最大值的位置可以填(le h+1)的数。
可以发现(A,B)两矩阵的填入无交，那么分为两个阶段转移即可，预处理幂之后时间复杂度为(O(nmk(n+m)))。
#include<cstdio>
const int N=107;
int n,m,k,P,f[N][N][N],C[N][N],pw[N][N];
void inc(int&a,int b){a+=b-P,a+=a>>31&P;}
int sub(int a,int b){return a-=b,a+(a>>31&P);}
int mul(int a,int b){return 1ll*a*b%P;}
int main()
{
scanf("%d%d%d%d",&n,&m,&k,&P),f[0][0][0]=1;
for(int i=0;i<=100;++i) for(int j=C[i][0]=1;j<=i;++j) inc(C[i][j]=C[i-1][j],C[i-1][j-1]);
for(int i=0;i<=100;++i) for(int j=pw[i][0]=1;j<=100;++j) pw[i][j]=mul(pw[i][j-1],i);
for(int h=1;h<=k;++h)
{
for(int i=0;i<=n;++i)
for(int j=0;j<=m;++j)
for(int t=mul(pw[k-h+1][j],sub(pw[h][m-j],pw[h-1][m-j])),s=1,l=i;l<=n;++l)
inc(f[h][l][j],mul(mul(s,C[n-i][l-i]),f[h-1][i][j])),s=mul(s,t);
for(int i=0;i<=n;++i)
for(int j=m;~j;--j)
for(int t=mul(pw[h][n-i],sub(pw[k-h+1][i],pw[k-h][i])),s=1,l=j+1;l<=m;++l)
s=mul(s,t),inc(f[h][i][l],mul(mul(s,C[m-j][l-j]),f[h][i][j]));
}
printf("%d",f[k][n][m]);
}