题目传送门：P10004 [集训队互测 2023] Permutation Counting 2。
考虑二项式反演，设 \(f_{i,j}\) 表示钦定 \(i\) 个位置为原排列上升位，\(j\) 个位置为逆排列上升位的方案数。
可以发现我们相当于钦定了原排列有 \(n-i\) 个上升连续段，逆排列有 \(n-j\) 个上升连续段，进一步观察发现，显然对于一个原排列的上升连续段的元素位置一定也上升，相当于一个原排列的连续段划分成了若干个逆排列的连续段，设 \(a_{i,j}\) 表示原排列第 \(i\) 个连续段有 \(a_{i,j}\) 个在逆排列的第 \(j\) 段，显然 \(a\) 需满足每一行每一列都不为 \(0\) 且总和为 \(n\)，对这个东西计数。
设 \(g_{i,j}\) 为原排列有 \(i\) 个上升连续段，逆排列有 \(j\) 个上升连续段的方案数，再次考虑容斥，枚举行列钦定总和不为 \(0\) 的位置，然后剩下的位置的计数就是一个插板法。
那么 \(g_{i,j} = \sum\limits_{x=0}^i (-1)^{i-x} {i\choose x} \sum\limits_{y=0}^j (-1)^{j-y} {j\choose y} {n+xy-1 \choose xy-1}\)，\(f_{i,j}=g_{n-i,n-j}\)。
最后的答案直接二项式反演即可 \(ans_{i,j}=\sum\limits_{x=i}^n (-1)^{x-i} {x\choose i} \sum\limits_{y=j}^n (-1)^{y-j} {y\choose j} f_{x,y}\)，直接做是 \(O(n^4)\) 的，考虑优化。
设 \(s_{x,j}= \sum\limits_{y=0}^j (-1)^{j-y} {j\choose y} {n+xy-1\choose xy-1}\)，那么 \(g_{i,j}=\sum\limits_{x=0}^i (-1)^{i-x} {i\choose x} s_{x,j}\)。
设 \(h_{x,j}=\sum\limits_{y=j}^n (-1)^{y-j} {y\choose j} f_{x,y}\)，那么 \(ans_{i,j}=\sum\limits_{x=i}^n (-1)^{x-i} {x\choose i} h_{x,j}\)。
这样就是 \(O(n^3)\) 的，有点卡常，可以预处理组合数减少取模量。
#include<bits/stdc++.h>
#define int long long
#define double long double
using namespace std;
const int N=510;
inline int read(){char c=getchar();int f=1,ans=0;while(c<48||c>57) f=(c==45?f=-1:1),c=getchar();while(c>=48&&c<=57) ans=(ans<<1)+(ans<<3)+(c^48),c=getchar();return ans*f;
}
int n,mod;
int fac[N*N],invfac[N*N],f[N][N],pw[N*N],s[N][N],h[N][N],c[N][N],t[N][N];
inline void add(int &x,int y){x+=y;if (x>=mod) x-=mod;}
void exgcd(int a,int b,int &x,int &y){if (b==0) x=1,y=0;else exgcd(b,a%b,y,x),y-=a/b*x;}
inline int inv(int a){int x,y;exgcd(a,mod,x,y);return (x%mod+mod)%mod;}
inline int C(int n,int m){return fac[n]*invfac[m]%mod*invfac[n-m]%mod;}
main(){n=read(),mod=read();fac[0]=invfac[0]=pw[0]=1;for (int i=1;i<=n*n+n;i++) fac[i]=fac[i-1]*i%mod,invfac[i]=invfac[i-1]*inv(i)%mod,pw[i]=pw[i-1]*(mod-1)%mod;c[0][0]=1;for (int i=1;i<=n;i++){c[i][0]=1;for (int j=1;j<=i;j++) c[i][j]=(c[i-1][j]+c[i-1][j-1])%mod;}for (int x=1;x<=n;x++) for (int y=1;y<=n;y++) t[x][y]=C(n+x*y-1,x*y-1);for (int x=0;x<=n;x++) for (int j=1;j<=n;j++) for (int y=0;y<=j;y++) add(s[x][j],pw[j-y]*c[j][y]%mod*t[x][y]%mod);for (int i=1;i<=n;i++) for (int j=1;j<=n;j++) for (int x=0;x<=i;x++) add(f[i][j],pw[i-x]*c[i][x]%mod*s[x][j]%mod);for (int x=0;x<n;x++) for (int j=0;j<n;j++) for (int y=j;y<n;y++) add(h[x][j],pw[y-j]*c[y][j]%mod*f[n-x][n-y]%mod);for (int i=0;i<n;i++,puts("")) for (int j=0;j<n;j++){int ans=0;for (int x=i;x<n;x++) add(ans,pw[x-i]*c[x][i]%mod*h[x][j]%mod);printf("%lld ",ans);}return 0;
}