给定正整数 \(n,k\) 和质数 \(p\)，定义长为 \(k\) 的正整数序列 \(\{a_1,\cdots,a_k\}\) 是好的 \(\iff\forall i\in[1,k],a_i\le i\)。
对每个 \(\text{pos}\in[1,nk],\text{val}\in[1,k]\)，求使得 \(b_{\text{pos}}=\text{val}\) 且可以划分为 \(n\) 个长为 \(k\) 的好子序列的长为 \(nk\) 的正整数序列 \(\{b_1,\cdots,b_{nk}\}\) 的数量 \(\bmod p\)。
\(n,k\le 20\)，\(10^8<p<10^9\)。
考虑判断一个序列是不是 amazing 的，直接从后往前贪心放，放到能放的子序列里面位置最前的。
设每个子序列剩余位置为 \(c_1,\cdots,c_n\)，初始时全为 \(k\)，贪心过程即找到不小于 \(b_i\) 的第一个位置然后减 \(1\)。
考虑进行一个转置：将剩余位置排成一个矩形，维护左边界，就转化为 \(k\) 个不同的球在数轴上，初始位置全为 \(0\)，每次操作选择一个位置不在 \(n\) 的球，将其移到下一个位置。
\(b_{\text{pos}}=\text{val}\) 的限制即为第 \(nk-\text{pos}+1\) 次操作将位置为 \(x\) 的球放到 \(x+1\)，其中 \(\text{pre}_x<\text{val}\le\text{pre}_{x+1}\)，\(\text{pre}_x\) 表示位置 \(0,1,\cdots,x-1\) 上的球的个数。
设这次操作之前第 \(i\) 个球的位置为 \(p_i\)，则 \(p_1+\cdots+p_k=\text{pos}-1\)，方案数为 \(\dfrac{(i-1)!(nk-i)!(n-p_?)}{\prod p_j!(n-p_j)!}\)，其中 \(p_?\) 表示这次操作的球的位置。
然后就可以 dp 了，先枚举 \(\text{val}\)，设 \(\text{cnt}_i\) 表示有多少个 \(p_j\) 为 \(i\)，从小到大枚举 \(i\)，将 \(\text{pre}_i=\sum\limits_{j<i}\text{cnt}_j\) 和 \(\sum\limits_{j<i}j\cdot\text{cnt}_j\) 记到状态里，就可以知道 \(p_?\) 是哪一个了，用 dp 数组即可算出每个 \(\text{pos}\) 的答案。时间复杂度 \(O(n^2k^4)\)。
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 21, M = N*N;
int n, k, mod, pw[N][N], fac[M], inv[M], ans[M][N], f[N][M], g[N][M];
int ksm(int a, int b){
int res = 1;
for(;b;b >>= 1, a = (LL)a * a % mod)
if(b & 1) res = (LL)res * a % mod;
return res;
}
int main(){
ios::sync_with_stdio(false);
cin >> n >> k >> mod; *fac = 1;
for(int i = 1;i < M;++ i) fac[i] = (LL)fac[i-1] * i % mod;
inv[M-1] = ksm(fac[M-1], mod-2);
for(int i = M-1;i;-- i) inv[i-1] = (LL)inv[i] * i % mod;
for(int i = 0;i < N;++ i)
for(int j = pw[i][0] = 1;j < N;++ j)
pw[i][j] = (LL)pw[i][j-1] * inv[i] % mod;
for(int val = 1;val <= k;++ val){
memset(f, 0, sizeof f); **f = 1;
for(int i = 0;i < n;++ i){
memset(g, 0, sizeof g);
for(int j = 0;j <= k;++ j)
for(int l = 0;l <= j*(i-1);++ l) if(f[j][l])
for(int t = 0;t <= k-j;++ t)
if(j < val && val <= j+t) g[j+t][l+i*t] = (g[j+t][l+i*t] + LL(n-i) * f[j][l] % mod * inv[t] % mod * pw[i][t] % mod * pw[n-i][t]) % mod;
else g[j+t][l+i*t] = (g[j+t][l+i*t] + (LL)f[j][l] * inv[t] % mod * pw[i][t] % mod * pw[n-i][t]) % mod;
memcpy(f, g, sizeof f);
}
for(int i = 1;i <= n*k;++ i)
for(int j = val;j <= k;++ j)
if(i-1 >= (k-j)*n) ans[n*k-i+1][val] = (ans[n*k-i+1][val] + (LL)f[j][i-1-(k-j)*n] * inv[k-j] % mod * pw[n][k-j]) % mod;
}
for(int i = 1;i <= n*k;++ i)
for(int j = 1;j <= k;++ j)
printf("%lld%c", (LL)ans[i][j]*fac[k]%mod*fac[i-1]%mod*fac[n*k-i]%mod, " \n"[j==k]);
}