一道需要一些猜结论技巧的中档题。
首先突破口在于排列长度恰好等于不是额外输入的某个数 \(k\) 而是 LDS 与 LIS 的乘积，这显然启示我们去找一些性质。根据 dilworth 定理，最长反链等于最小链覆盖，故 LIS 的长度，就是最少需要的递减数列的个数使得每个元素被覆盖至少一次，而每个递减数列长度不超过 \(m\)，因此 \(n\) 个递减数列长度不超过 \(nm\)。换句话，\(nm\) 是 \(\text{LIS}=n,\text{LDS}=m\) 的排列长度的最小值。
先考虑怎么计算 \(\text{LIS}=n,\text{LDS}=m,\text{len}=nm\) 的排列个数。大胆猜一个结论：对于符合要求的排列，将其划分为 \(n\) 个递减数列的方案是唯一的。并且如果我们假设 \(a_{i,j}\) 为第 \(i\) 个递减数列中第 \(j\) 个元素在排列中的位置，\(p_{i,j}\) 为第 \(i\) 个递减数列中第 \(j\) 个元素的值，那么有：
- \(a_{i,j}\le a_{i,j+1},a_{i,j}\le a_{i+1,j}\)
- \(p_{i,j}\ge p_{i,j+1},p_{i,j}\le p_{i+1,j}\)。
发现 \(a,p\) 的安排是两个完全不相干的子问题，因此分别求解方案数然后乘起来就是排列总数。而两者方案数显然是完全相同的，写个阶乘的暴力验证一下也会发现答案确实是完全平方数。因此只用知道怎么计算符合条件的 \(a\) 的个数即可。而这个过程显然相当于每次新填一个数的时候，已经填的数和没填的数的分界线是一条从左下角到右上角，且只能向右或向上的折线。这样的折线最多只有 \(\dbinom{n+m}{n}\) 种，当 \(nm\le 100\) 时该值大概最多是 \(2\times 10^5\) 级别的，因此直接写个状压 dp 记录折线信息即可。
知道怎么算总数以后，计算 \(f(pos,val)\) 的部分是 trivial 的——枚举这个位置上的值是第 \(x\) 个下降序列的第 \(y\) 个元素，那么对 \(a,p\) 的限制相当于强制要求 \(a_{x,y}=v\)，预处理 \(F_{x,y,v}\) 表示强制令 \(a_{x,y}=v\) 的符合要求的矩阵的方案数，这个可以直接枚举左边部分的折线然后 DP 算一下左右部分的方案数乘起来得到。
const int MAXN=100;
const int MAXC=184756;
int n,m,mod,rv,f[MAXN+5][MAXN+5][MAXN+5];
int calc(int x,int y,int num){
if(rv)return f[y][x][num];
else return f[x][y][num];
}
bool check(vector<int>v){
for(int i=1;i<v.size();i++)if(v[i]>v[i-1])return 0;
return 1;
}
map<vector<int>,int>id;
vector<int>cur,v[MAXC+5];
int idcnt=0,dp[MAXC+5];
void dfs(int x,int lst){
if(x==n+1){
reverse(cur.begin(),cur.end());
id[cur]=++idcnt;v[idcnt]=cur;
reverse(cur.begin(),cur.end());return;
}
for(int i=lst;i<=m;i++)cur.pb(i),dfs(x+1,i),cur.ppb();
}
int main(){
scanf("%d%d%d",&n,&m,&mod);
if(n>m)swap(n,m),rv=1;
dfs(1,0);dp[1]=1;
for(int i=1;i<=idcnt;i++){
vector<int>cur=v[i];
for(int j=0;j<n;j++)if(cur[j]<m){
cur[j]++;
if(check(cur)){
int nid=id[cur];
dp[nid]=(dp[nid]+dp[i])%mod;
}
cur[j]--;
}
}
for(int i=1;i<=idcnt;i++){
vector<int>cur=v[i];
for(int j=0;j<n;j++)if(cur[j]<m){
cur[j]++;
if(check(cur)){
int nid=id[cur],sum=0;for(int k=0;k<n;k++)sum+=cur[k];
vector<int>rst(n);
for(int k=0;k<n;k++)rst[k]=m-cur[n-1-k];
f[j+1][cur[j]][sum]=(f[j+1][cur[j]][sum]+1ll*dp[i]*dp[id[rst]])%mod;
}cur[j]--;
}
}
if(rv)swap(n,m);
for(int i=1;i<=n*m;i++)for(int j=1;j<=n*m;j++){
int res=0;
for(int x=1;x<=n;x++)for(int y=1;y<=m;y++)
res=(res+1ll*calc(x,y,i)*calc(x,m-y+1,j))%mod;
printf("%d%c",res," \n"[j==n*m]);
}
return 0;
}