把排列写成一条路径 \(p_1\to p_2\to\cdots\to p_n\)。那么 \([p_i<p_{i+1}]\) 就是第 \(i\) 步往右走，\([p^{-1}_i<p^{-1}_{i+1}]\) 是 \(i\) 要先于 \(i+1\) 访问。
如果我们已知了 \(p^{-1}_{i}<p^{-1}_{i+1}\)，此时路径形如 \(\cdots\to i\to\cdots\to i+1\to\cdots\)。可以发现，此时我们可以把 \(i\) 和 \(i+1\) 缩在一起，并允许它走两遍。第一次选择 \(i\)，第二次选择 \(i+1\)。
考虑钦定若干组 \(p^{-1}_{i}<p^{-1}_{i+1}\)，此时序列会被缩成 \(k\) 段，长度依次为 \(c_1,c_2,\cdots,c_k\)。因为每段内必须从小往大走，所以其实只需要考虑走过的点所属连续段。
不妨设走过的点所属连续段为 \(v_1,v_2,\cdots,v_n\)。容易发现，此时 \(\sum[p_i<p_{i+1}]\) 相当于 \(\sum[v_i\le v_{i+1}]\)。那么我们想要解决的就是把 \(c_i\) 个 \(i\) 重排，数 \(\sum[v_i\le v_{i+1}]=x\) 的方案数。这个可以 \(O(n^5)\) dp。
$O(n^5)$ 的实现
#include <bits/stdc++.h>
using namespace std;
using ll = long long;int mod;
struct Barrett {__uint128_t t;Barrett(int mod = 2) {t = ((__uint128_t)1 << 64) / mod;}ll operator () (ll x) {x -= ((t * x) >> 64) * mod;return x -= (x >= mod) * mod;}
}MOD;
void Add(int &x, ll y) { x = MOD(x + y); }const int kN = 505;
int n;
int C[kN][kN], mul[kN];
int f[kN][kN][kN], g[kN][kN][kN];
int ans[kN][kN];void InitComb(int N = kN - 1) {mul[0] = 1;for(int i = 1; i <= N; i++) {mul[i] = MOD((ll)mul[i - 1] * i);}for(int i = 0; i <= N; i++) {C[i][0] = 1;for(int j = 1; j <= i; j++) {C[i][j] = MOD(C[i - 1][j] + C[i - 1][j - 1]);}}
}void PreWork() {for(int i = 1; i <= n; i++) {for(int j = 1; i + j <= n + 1; j++) {for(int k = 0; k <= min(i, j); k++) {for(int c = 0; c <= k; c++) {int val = MOD((ll)C[k][c] * C[j - c + i - 1][i]);Add(g[i][j][k], (c & 1) ? mod - val : val);}}}}
}int main() {// freopen("1.in", "r", stdin);// freopen("1.out", "w", stdout);ios::sync_with_stdio(0), cin.tie(0);cin >> n >> mod;MOD = Barrett(mod);InitComb();PreWork();for(int i = 1; i <= n; i++) f[1][i][i - 1] = 1;for(int i = 1; i <= n; i++) {for(int s = 1; s <= n; s++) {for(int j = 0; j <= s; j++) {if(!f[i][s][j]) continue;ll val = f[i][s][j];for(int c = 1; s + c <= n; c++) {for(int t = 0; t <= min(c, j + 1); t++) {ll coe = MOD((ll)C[j + 1][t] * g[c][s - j + t][t]);Add(f[i + 1][s + c][j + c - t], val * coe);}}}}for(int c = 0; c < n; c++) Add(ans[n - i][c], f[i][n][c]);}for(int i = 0; i < n; i++, cout << "\n") {for(int j = 0; j < n; j++) {int sum = 0;for(int k = i; k < n; k++) {int val = MOD((ll)ans[k][j] * C[k][i]);Add(sum, ((k - i) & 1) ? mod - val : val);}cout << sum << " ";}}return 0;
}
实际上可以发现，这样直接 dp 状态数已经是三方了，转移看起来也不太能 \(O(1)\)。
对原路径 \(p_1\to\cdots\to p_n\) 分析。此时 \(k\) 段就是 \(k\) 个区间 \([l_i,r_i]\)。
考虑从前往后扫，维护相对大小关系。一个想法是 \(f_{i,j,k}\) 表示前 \(i\) 个位置，\(j\) 个区间，钦定 \(k\) 个上升的方案数，每次加入区间。但是这样转移需要枚举新增的递增长度和区间数量，似乎又只能五方。
另一个想法是，枚举区间数 \(c\)。然后 \(f_{i,j}\) 表示前 \(i\) 个位置，钦定 \(j\) 个上升。这样可以总复杂度 \(O(n^4)\) 地 dp 出来。然后容斥一下空区间和把一个区间拆成多个，这样总复杂度也是 \(O(n^4)\)。
$O(n^4)$ 的实现
#include <bits/stdc++.h>
using namespace std;
using ll = long long;int mod;
struct Barrett {__uint128_t t;Barrett(int mod = 2) {t = ((__uint128_t)1 << 64) / mod;}ll operator () (ll x) {x -= ((t * x) >> 64) * mod;return x -= (x >= mod) * mod;}
}MOD;
void Add(int &x, ll y) { x = MOD(x + y); }const int kN = 1005;
int n;
int C[kN][kN], mul[kN];
int f[kN][kN];
int ans[kN][kN];void InitComb(int N = kN - 1) {mul[0] = 1;for(int i = 1; i <= N; i++) {mul[i] = MOD((ll)mul[i - 1] * i);}for(int i = 0; i <= N; i++) {C[i][0] = 1;for(int j = 1; j <= i; j++) {C[i][j] = MOD(C[i - 1][j] + C[i - 1][j - 1]);}}
}int main() {// freopen("1.in", "r", stdin);// freopen("1.out", "w", stdout);ios::sync_with_stdio(0), cin.tie(0);cin >> n >> mod;MOD = Barrett(mod);InitComb();for(int c = 1; c <= n; c++) {memset(f, 0, sizeof(f));f[0][0] = 1;for(int i = 0; i < n; i++) {for(int j = 0; j <= i; j++) {ll val = f[i][j];for(int k = 1; k + i <= n; k++) {Add(f[i + k][j + k - 1], val * C[c + k - 1][k]);}}}memcpy(ans[c], f[n], sizeof(ans[c]));}for(int i = n; i; i--) {for(int j = 0; j < n; j++) {int sum = 0;for(int k = i; k; k--) {ll val = MOD((ll)ans[k][j] * C[i][k]);Add(sum, ((i - k) & 1) ? mod - val : val);}ans[i][j] = sum;}}for(int i = n; i; i--) {for(int j = 0; j < n; j++) {int sum = 0;for(int k = i; k; k--) {ll val = MOD((ll)ans[k][j] * C[n - k][i - k]);Add(sum, ((i - k) & 1) ? mod - val : val);}ans[i][j] = sum;}}for(int i = 1; i <= n; i++) {for(int j = 0; j < n; j++) {int sum = 0;for(int k = j; k < n; k++) {ll val = MOD((ll)ans[i][k] * C[k][j]);Add(sum, ((k - j) & 1) ? mod - val : val);}ans[i][j] = sum;}}for(int i = 0; i < n; i++, cout << "\n") {for(int j = 0; j < n; j++) {cout << ans[n - i][j] << " ";}}return 0;
}
之后就是简单的了。这个做法的瓶颈在于对 \(f\) 的 dp。转移是 \(f_{i,j}\times\dbinom{c+k-1}{k}\to f_{i+k,j+k-1}\)。
令 \(g_{i,i-j}=f_{i,j}\)，那么转移可以改写成 \(g{i,i-j}\times\dbinom{c+k-1}{k}\to g_{i+k,i-j+1}\)。此时就是卷积的形式。写出转移的 OGF 为 \(F(x)=\sum\limits_{n\ge 1}x^n\dbinom{c+n-1}{n}=\dfrac{1}{(1-x)^c}-1\)。
那么我们要算形如 \([x^m]F^n(x)\)，拆一下括号即可。总复杂度 \(O(n^3)\)。
#include <bits/stdc++.h>
using namespace std;
using ll = long long;int mod;
struct Barrett {__uint128_t t;Barrett(int mod = 2) {t = ((__uint128_t)1 << 64) / mod;}ll operator () (ll x) {x -= ((t * x) >> 64) * mod;return x -= (x >= mod) * mod;}
}MOD;
void Add(int &x, ll y) { x = MOD(x + y); }
int Pow(int x, int y) {int b = x, r = 1;for(; y; y /= 2, b = MOD((ll)b * b)) {if(y & 1) r = MOD((ll)r * b);}return r;
}const int kN = 505, kL = 3e5 + 5;
int n;
int mul[kL], imul[kL], cm[kL];
int C[kN][kN];
int ans[kN][kN];void InitComb(int N = kL - 1, int M = kN - 1) {mul[0] = 1;for(int i = 1; i <= N; i++) {mul[i] = MOD((ll)mul[i - 1] * i);}imul[N] = Pow(mul[N], mod - 2);for(int i = N - 1; ~i; i--) {imul[i] = MOD((ll)imul[i + 1] * (i + 1));}for(int i = 0; i <= M; i++) {C[i][0] = 1;for(int j = 1; j <= i; j++) {C[i][j] = MOD(C[i - 1][j] + C[i - 1][j - 1]);}}
}
int Comb(int n, int m) {if((n < m) || (m < 0)) return 0;return MOD(MOD((ll)mul[n] * imul[m]) * imul[n - m]);
}int main() {// freopen("1.in", "r", stdin);// freopen("1.out", "w", stdout);ios::sync_with_stdio(0), cin.tie(0);cin >> n >> mod;MOD = Barrett(mod);InitComb();for(int i = 0; i < kL; i++) cm[i] = Comb(i, n);for(int c = 1; c <= n; c++) {for(int i = 0; i < n; i++) {for(int j = 0; j <= n - i; j++) {ll val = MOD((ll)C[n - i][j] * cm[c * (n - i - j) + n - 1]);Add(ans[c][i], (j & 1) ? mod - val : val);}}}for(int i = n; i; i--) {for(int j = 0; j < n; j++) {int sum = 0;for(int k = i; k; k--) {ll val = MOD((ll)ans[k][j] * C[i][k]);Add(sum, ((i - k) & 1) ? mod - val : val);}ans[i][j] = sum;}}for(int i = n; i; i--) {for(int j = 0; j < n; j++) {int sum = 0;for(int k = i; k; k--) {ll val = MOD((ll)ans[k][j] * C[n - k][i - k]);Add(sum, ((i - k) & 1) ? mod - val : val);}ans[i][j] = sum;}}for(int i = 1; i <= n; i++) {for(int j = 0; j < n; j++) {int sum = 0;for(int k = j; k < n; k++) {ll val = MOD((ll)ans[i][k] * C[k][j]);Add(sum, ((k - j) & 1) ? mod - val : val);}ans[i][j] = sum;}}for(int i = 0; i < n; i++, cout << "\n") {for(int j = 0; j < n; j++) {cout << ans[n - i][j] << " ";}}return 0;
}
)
