給定正整數 \(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\)

設這次操作之前第 \(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\)

#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]);
}