对于一个序列 $a$，似乎很难计算其移除标记的方案数 $f(a)$。
不妨换一个思路，对于一个移除标记的方案，去计算有多少序列能够产生这种方案。

具体的说，记 $b_i$ 为第 $i$ 次操作中移出的 token 位置（如果 $a_i=0$ 没有移出 token 则记 $b_i=0$），那么有限制所有 $>0$ 的 $b_i$ 互不相同。
此时如果固定了数组 $b$，那这个 $b$ 对应的序列 $a$ 的数量即为 $\prod\limits_{b_i>0}b_i$，因为每个 $a_i(a_i\ne0)$ 的取值都是 $[1,b_i]$。

不妨把问题放到棋盘上。棋盘上有 $n\times n$ 个点，要求在其中放不超过 $n$ 个车，要求车不能互相攻击，并且每个车所在的位置 $(x,y)$ 必须满足 $x\ge y$。
记一种放置方法的权值为所有车所在 $y$ 坐标之积，求所有放置方法的权值之和。

如果还是按照之前的想法，第 $i$ 次 dp 更新 $x=i$ 的车，仍然非常困难，因为无法记录哪些 $y$ 被用过。
不妨参考一下 ferrers 图，转一下棋盘，第 $i$ 次 dp 更新 $y=i$ 的车。
但是还是很困难，因为每次能放的 $x$ 会越来越小，无法记录之前有多少个车的 $x$ 对后面有影响。
不妨倒着做，第 $i$ 次 dp 更新 $y=i$ 的车，但是转移是从 $n$ 到 $1$。
记录 $f_{i,j}$ 表示已经决定 $y\in[i,n]$ 的车的放置情况，并且一共放了 $j$ 个车。
转移有两种情况：

  • 不放车，权值继承 $f_{i+1,j}$；
  • 放车，由于 $x\ge y=i$，因此 $x$ 有 $n-i+1$ 种选法，由于已经放了 $j-1$ 个车，因此还有 $(n-i+1)-(j-1)$ 种，权值的贡献就是 $f_{i+1,j-1}\times((n-i+1)-(j-1))\times i$。

最终答案就是 $\sum\limits_{j=0}^nf_{1,j}$。

#include<bits/stdc++.h>
#define siz(x) int((x).size())
#define all(x) std::begin(x),std::end(x)
#define fi first
#define se second
using namespace std;
using unt=unsigned;
using loli=long long;
using lolu=unsigned long long;
using pii=pair<int,int>;
mt19937_64 rng(random_device{}());
constexpr int N=5007;
int P;
struct mint{
	int d;
	mint()=default;
	mint(int x):d(x){}
	friend std::istream&operator>>(std::istream&x,mint&y){return x>>y.d;}
	friend std::ostream&operator<<(std::ostream&x,mint y){return x<<y.d;}
	friend mint operator+(mint x,mint y){return (x.d+=y.d)<P?x.d:x.d-P;}
	mint&operator+=(mint z){return (d+=z.d)<P?d:d-=P,*this;}
	friend mint operator-(mint x,mint y){return (x.d-=y.d)<0?x.d+P:x.d;}
	mint&operator-=(mint z){return (d-=z.d)<0?d+=P:d,*this;}
	friend mint operator*(mint x,mint y){return int(1ll*x.d*y.d%P);}
	mint&operator*=(mint z){return d=int(1ll*d*z.d%P),*this;}
	static mint qpow(int x,int y=P-2){int z=1;for(;y;y>>=1,x=int(1ll*x*x%P))if(y&1)z=int(1ll*x*z%P);return z;}
	friend mint operator/(mint x,mint y){return x*=qpow(y.d);}
	mint&operator/=(mint z){return (*this)*=qpow(z.d);}
	friend mint operator^(mint x,mint y){return qpow(x.d,y.d);}
	mint&operator^=(mint z){return *this=qpow(d,z.d);}
	mint operator()(mint z)const{return qpow(d,z.d);}
	mint&operator[](mint z){return *this=qpow(d,z.d);}
	mint inv()const{return qpow(d);}
	mint pow(mint z)const{return qpow(d,z.d);}
	int operator+()const{return d;}
	mint operator-()const{return d?P-d:0;}
	int operator~()const{return ~d;}
};
mint operator""_m(lolu x){return mint(int(x%P));}
int n;
mint f[N][N];
signed main(){
//	freopen(".in","r",stdin);
//	freopen(".out","w",stdout);
	ios::sync_with_stdio(false);cin.tie(nullptr);
	int T;cin>>T;while(T--){
		cin>>n>>P;
		for(int i=1;i<=n+1;i++)
			memset(f[i],0,sizeof(int)*(n+2));
		f[n+1][0]=1;
		for(int i=n;i>=1;i--)for(int j=0;j<=n-i+1;j++){
			f[i][j]=f[i+1][j];
			if(j)f[i][j]+=f[i+1][j-1]*((n-i+1)-(j-1))*i;
		}
		mint ans=0;
		for(int j=0;j<=n;j++)ans+=f[1][j];
		cout<<ans<<'\n';
	}
	return 0;
}