很神奇的 \(dp\) 题。
我们观察到 \(k^2\) = \(\binom{k}{2}×2\) + \(k\)。
实际上我们可以把问题转化为求一个数 \(x\) 有几次出现以及在每个序列任选两个 \(x\) 的总方案数。

定义：

\(f_{i,j}\) ： 序列的前缀最大值为 \(j\)，长度为 \(i\) 的方案数。
\(g_{i,j}\) ： 在一个前缀最大值为 \(j\) 的序列后再选 \(i\) 个数的方案数。

转移：

• \(f_{i,j}=f_{i-1,j}×j+f_{i-1,j-1}\)
  这个应该很好理解，就是选数。
• \(g_{i,j}=g_{i-1,j}×j+g_{i-1,j+1}\)
  前半部分实际上是在该序列再加一个小于等于曾经的前缀最大值的数，后半部分的实际意义是在曾经的前缀最大值比自己大1的 \(g\) 数组的一个长度为 \(i-1\) 的序列的前端加上一个 \(j\) 。可以发现，转移不重不漏。

答案更新

对于一个在 \(i\) 位置上的 \(x\)，其对答案的贡献是
\(\sum_{y=x}^{n}f_{i-1,y}（g_{n-i,y}+2(n-i)g_{n-i-1,y}）+f_{i-1,x-1}(g_{n-i,x}+2(n-i)g_{n-i-1,x})\)
转移总体分为更新最大值和不更新最大值。
对于每部分的前一半是在找该位上的 \(x\) 共出现了几次，后半部分是从长度为 \(n-i\) 的序列后再选一个 \(x\) 的方案数，以此来找每个序列任选两个 \(x\) 的总方案数。转移时记录前缀和就行了。

code

#include <bits/stdc++.h>
#define re register
#define int long long
#define db double
#define pir make_pair
using namespace std; 
const int maxn=3010;
const int INF=1e18;
inline int read() {
   int s=0,w=1; char ch=getchar();
   while(ch<'0'||ch>'9') { if(ch=='-') w=-1;ch=getchar(); }
   while(ch>='0'&&ch<='9') { s=s*10+ch-'0'; ch=getchar(); }
   return s*w;
}
int n,mol,f[maxn][maxn],sum[maxn][maxn],g[maxn][maxn];
inline void add(int &x) { x= x>=mol? x-mol:x; }
signed main(void) {
   //freopen("cs.in","r",stdin);
   n=read(),mol=read();
   f[0][0]=1;
   for(re int i=1;i<=n;i++) for(re int j=1;j<=n;j++) f[i][j]=(f[i-1][j]*j%mol+f[i-1][j-1])%mol;
   for(re int i=1;i<=n;i++) g[0][i]=1;
   for(re int i=1;i<=n;i++) for(re int j=1;j<=n;j++) g[i][j]=(g[i-1][j]*j%mol+g[i-1][j+1])%mol;
   for(re int i=1;i<n;i++) {
   	for(re int y=n;y>=1;y--) {
   		sum[i][y]=(sum[i][y+1]+f[i-1][y]*(g[n-i][y]+2*(n-i)%mol*g[n-i-1][y]%mol)%mol)%mol;
   	}
   }
   for(re int y=n;y>=1;y--) 
   	sum[n][y]=(sum[n][y+1]+f[n-1][y]*g[n-n][y]%mol)%mol;
   for(re int x=1,res;x<=n;x++) {
   	res=0;
   	for(re int i=1;i<n;i++) {
   		add(res+=sum[i][x]);
   		add(res+=(f[i-1][x-1]*(g[n-i][x]+2*(n-i)%mol*g[n-i-1][x]%mol)%mol));
   	}
   	add(res+=(sum[n][x]+f[n-1][x-1]*g[n-n][x]%mol)%mol);
   	printf("%lld ",res);
   }
   return 0;
}