Problem Link

首先考虑刻画等价类关系，限定两个子树同构较简单，而限定两个子树不同构则有一定难度。

那么我们进行集合不相等容斥，我们先按子树同构将每个点划分等价类，然后容斥钦定若干对等价类同构。

此时我们会合并这些被钦定边的等价类，假设我们把原来的 $k$ 个等价类合并成一个，则系数为所有 $k$ 个点的连通图 $-1$ 的边数次方之和，这是经典问题，结果是 $(-1)^{k-1}(k-1)!$。

那么我们考虑枚举容斥后的每个等价类大小，设大小为 $i$ 的等价类有 $a_i$ 个，则分配标号的方案数为 $\dfrac{(n-1)!}{\prod a_i!\times (i!)^{a_i}}$，其中 $(n-1)!$ 是因为要钦定大小为 $1$ 的点为根。

然后我们考虑如何连接树结构，显然每个等价类 $u$ 中点的父亲都处于更小的等价类 $v$，且如果存在 $fa_u\to v$，则 $v$ 的每个点都要从一个 $u$ 中点连接。

也就是每个 $u$ 的父亲应该是若干个等价类 $v$，每个等价类中的每个点恰好是 $u$ 中一个点的父亲，即 $\sum siz_v=siz_u$。

那么我们考虑如何在已知 $a_1\sim a_{k-1}$ 的连边情况时确定 $a_k$ 个大小为 $k$ 的等价类连边方案。

首先我们把这些连通块分成两类，第一类是树根，这些点的父是大小更小等价类，第二类是非树根，这些点父亲是一些其他的大小为 $k$ 的等价类。

先计算生成一个一类连通块的方案。

首先考虑单个连通块的 EGF $A(x)=\prod_{i<k}\left(\sum_j\left(\dfrac{x^{j}}{j!}\right)^i\right)^{a_i}$，其中 $j$ 是枚举向该连通块的边数。

加上 $q$ 的权值得到 $B(x)=q(A(x)-1)$，因为要调整常数项，然后我们要按照最开始说的集合不相等容斥，最终答案就是 $[i!x^i]\sum_{k>1}\dfrac{(-1)^{k-1}(k-1)!B^k(x)}{k!}=[i!x^i]\ln(1+B(x))$，

其中 $i!$ 的贡献和计算标号时的 $i!^{-a_i}$ 可以抵消。

而我们还要计算 $c$ 个点中有 $k$ 个根，生成森林的方案数，即 $k\times c^{c-k-1}$。

那么这样可以用 $\mathcal O(\pi(n)\mathrm{poly}(n))$ 的复杂度实现，其中 $\pi$ 是分拆数。

考虑如何优化 $\pi (n)$ 级别的状态量。

我们尝试对这些状态进行一定的 dp，然后把一些较简单的过程直接处理掉。

首先可以解决的就是所有叶子节点，对于是叶子的等价类，一定不会对更大的等价类连边产生影响，所以我们只要记录非叶子节点的等价类结构。

更具体的，记录向更大等价类有连边的点集，很显然这样的点数 $\le \dfrac n2$，状态数变为 $\pi(\frac n2)$ 级别，状态中还要记录当前已经加入的总点数。

实现时我们需要按等价类大小逐层加点，加入一层点之后需要删掉一些不再向后连边的点，为了保证我们的删点过程是正确的，还要再进行一轮相同操作，删掉若干个点并乘以 $-1$ 的系数使贡献正确。

时间复杂度 $\mathcal O(\pi(\frac n2)\mathrm{poly}(n))$。

考虑更优秀的算法，还要进一步压缩拆分数的状态常数。

我们考虑设定阈值 $B$，并且把所有点按子树大小是否 $>B$ 分类。

我们把这两种等价类分开处理，对于一个子树 $>B$ 的等价类，按照刚才的 dp 实现，对于子树 $\le B$ 的等价类，当 $B$ 很小时，爆搜出每一种有根树，并且枚举每种有根树是否作为子树出现过。

这里为了保证枚举的集合正确，也需要容斥，即枚举若干个有根树，并且钦定实际使用的有根树必须在这个集合中，可以得到若干种可用树的集合以及对应的系数（考虑上 $q$ 对这部分的贡献）。

然后对于子树大小 $>B$ 的点，分类成向后连边和未向后连边的点。

如果一个子树大小 $>B$ 的点向后有连边，那么按照刚才的方法记录进状态中，每次加边后删掉向后无边的点集并且容斥保证其正确性。

如果一个子树大小 $>B$ 的点向后无边，那么这个点要么在当前等价类中有儿子，否则必定由若干大小 $\le B$ 的子树合并成一棵大小 $>B$ 的子树。

那么在这部分中，我们实际上要考虑加入若干个大小为 $k$ 的等价类时，这些等价类所有 $\le B$ 的儿子。

这部分儿子不会影响 $q^s$ 计算，因为我们已经在外层枚举了可用有根树集合及其对应的系数，我们只要考虑其对加入总点数的贡献，计算出单个节点这样儿子关于大小的 OGF，预处理 $0\sim n$ 次幂即可。

此时对于原树的一些叶子节点，我们钦定其向其他 $>B$ 的子树无边，则这些点 $\le B$ 的儿子的大小之和 $\ge B$。

这依然不好计算，再进行一次容斥，枚举一些点 $\le B$ 的儿子大小之和 $<B$ 并向后无边，然后算出方案后再在剩下未被钦定的点中选一些钦定向后有边。

注意这里的组合系数是 $c$ 个点，选 $r$ 个点作为根并钦定 $c-p$ 个点作为叶子，我们枚举 $i$ 个被钦定叶子是根（单点），然后写出和式，可以用范德蒙德卷积优化掉 $i$ 的枚举。

此时被记录的点集子树大小 $\ge B+2$，则状态数是 $\pi(\frac{n}{B+2})$ 级别，取 $B=3$，此时 $\le B$ 的本质不同可用子树集合只有 $6$ 种，分别进行 dp 即可。

时间复杂度 $\mathcal O(\pi(\frac n5)\mathrm{poly}(n))$。

致谢：DeepSeek 老师的指导

代码：

#include<bits/stdc++.h>
#define ull unsigned long long
#define clr(a) memset(a,0,sizeof(a))
using namespace std;
namespace FastMod {
typedef __uint128_t uLL;
ull b,q,r; uLL m;
inline void init(const ull &B) { b=B,m=(uLL(1)<<64)/B; }
inline ull mod(const ull &a) {
	r=a-((m*a)>>64)*b;
	return r>=b?r-b:r;
}
}
int MOD;
struct mi {
	int x;
	mi(int X=0): x(X) {}
	#define o(x) FastMod::mod(x)
	mi operator +(const mi&o) { return x+o.x>=MOD?x+o.x-MOD:x+o.x; }
	mi operator -(const mi&o) { return x<o.x?x+MOD-o.x:x-o.x; }
	mi operator *(const mi&o) { return o(1ll*x*o.x); }
};
void operator +=(mi&x,const mi&o) { x=x+o; }
void operator -=(mi&x,const mi&o) { x=x-o; }
void operator *=(mi&x,const mi&o) { x=x*o; }
mi ksm(mi a,int b=MOD-2) { mi s=1; for(;b;a=a*a,b>>=1) if(b&1) s=s*a; return s; }
mi pw[105][105],pq[105],C[105][105],fac[105],ifac[105],ans[105],inv[105];
void Exp(mi*a,int n) {
	static mi b[105]; b[0]=1;
	for(int i=1;i<=n;++i) b[i]=0,a[i]*=i;
	for(int i=1;i<=n;++i) {
		for(int j=0;j<i;++j) b[i]+=b[j]*a[i-j];
		b[i]*=inv[i];
	}
	for(int i=0;i<=n;++i) a[i]=b[i];
}
void Ln(mi*a,int n) {
	static mi b[105];
	for(int i=0;i<=n;++i) b[i]=0;
	for(int i=1;i<=n;++i) {
		b[i]=a[i]*i;
		for(int j=1;j<i;++j) b[i]-=b[j]*a[i-j];
	}
	for(int i=0;i<=n;++i) a[i]=b[i]*inv[i];
}
void mul(mi*a,mi*b,int n) {
	static mi c[105];
	for(int i=0;i<=n;++i) {
		c[i]=0;
		for(int j=0;j<=i;++j) c[i]+=a[j]*b[i-j];
	}
	for(int i=0;i<=n;++i) a[i]=c[i];
}
const int MAXS=3005;
int n,bn,q,m,sz[MAXS],mx[MAXS],to[MAXS][25][25],ct[MAXS][25];
map <ull,int> id;
ull hw[MAXS][MAXS];
mt19937_64 rnd(time(0));
ull qhs(int*ar) {
	ull z=0;
	for(int i=0;i<=bn;++i) z^=hw[i][ar[i]];
	return z;
}
mi h[MAXS][105];
void init() {
	for(int i=0;i<=bn;++i) for(int j=0;j<=bn;++j) hw[i][j]=rnd();
	priority_queue<pair<int,vector<int>>>Q;
	Q.push({0,{}});
	while(Q.size()) {
		int vl=Q.top().first; auto u=Q.top().second;
		Q.pop(),sz[++m]=-vl; for(int x:u) ++ct[m][x];
		mx[m]=u.size()?u.back():0,id[qhs(ct[m])]=m,u.push_back(0);
		for(int&i=u.back()=max(1,mx[m]);i<=bn-sz[m];++i) Q.push({vl-i,u});
	}
	static mi p[105][105];
	for(int i=1;i<=bn;++i) {
		for(int j=0;j<=n/i;++j) p[i][j]=ksm(ifac[j],i);
		Ln(p[i],n/i);
	}
	static int a[25];
	for(int u=1;u<=m;++u) {
		for(int i=1;i<=bn;++i) {
			memcpy(a,ct[u],sizeof(a));
			for(int up=(bn-sz[u])/i+ct[u][i],&j=a[i]=0;j<=up;++j) to[u][i][j]=id[qhs(a)];
		}
		memcpy(a,ct[u],sizeof(a));
		for(int i=1;i<=bn;++i) if(a[i]) {
			for(int j=0;j<=n/i;++j) h[u][i*j]+=p[i][j]*a[i];
		}
		Exp(h[u],n),h[u][0]=1;
		for(int i=1;i<=n;++i) h[u][i]*=q;
		Ln(h[u],n);
	}
}
namespace lf {
int sz[4]={1,2,3,3},tr[4]={1,3,7,9};
mi w[4],z[16];
void init() {
	w[0]=w[1]=w[2]=1,w[3]=ksm(2);
	for(int s=1;s<16;++s) {
		int o=0,c=__builtin_popcount(s);
		for(int i=0;i<4;++i) if((tr[i]|s)==s) o|=1<<i;
		if(o) z[o]+=ksm(q,c)*ksm(1+MOD-q,4-c);
	}
}
}
mi f[MAXS][105],g[105][MAXS][105],w[105][105][105];
mi qL[105],qR[105],pL[105][105],pR[105][105];
signed main() {
	ios::sync_with_stdio(false),cin.tie(0),cout.tie(0);
	cin>>n>>q>>MOD,bn=n/5,FastMod::init(MOD);
	fac[0]=ifac[0]=pq[0]=1;
	for(int i=1;i<=n;++i) ifac[i]=ksm(fac[i]=fac[i-1]*i),pq[i]=pq[i-1]*q,inv[i]=ksm(i);
	for(int i=0;i<=n;++i) {
		C[i][0]=pw[i][0]=1;
		for(int j=1;j<=i;++j) C[i][j]=C[i-1][j]+C[i-1][j-1];
		for(int j=1;j<=n;++j) pw[i][j]=pw[i][j-1]*i;
	}
	lf::init(),init();
	for(int o:{1,3,7,11,15,9}) if(lf::z[o].x) {
		clr(f);
		for(int k=1,b;k<=(n-1)/4;++k) {
			b=n/k,clr(qR),qR[0]=1;
			for(int i:{0,1,2,3}) if(o>>i&1) {
				static mi t[105]; clr(t);
				for(int j=0,x=lf::sz[i];j<=b/x;++j) t[j*x]=ksm(lf::w[i],j*k)*ksm(ifac[j],k);
				mul(qR,t,b);
			}
			clr(qL),clr(pL),clr(pR),pL[0][0]=pR[0][0]=1;
			for(int i:{0,1,2}) qL[i]=qR[i];
			for(int i=1;i<=b;++i) {
				memcpy(pL[i],pL[i-1],sizeof(pL[i])),mul(pL[i],qL,b);
				memcpy(pR[i],pR[i-1],sizeof(pR[i])),mul(pR[i],qR,b);
			}
			clr(w);
			for(int r=1;r<=(k>1?b/4:1);++r) { //r: #root
				for(int c=r;c<=b;++c) for(int p=0;p<=c;++p) {
					//sum_i C(c,p)*C(c-p,r-i)*C(p,i)*p^{c-p-r+i}*p^{p-i-1}*i
					w[r][p][c-p]=mi((c-p)&1?MOD-1:1)*C[c][p]*C[c-1][r-1]*pw[p][c-r]*pq[c-r]*ifac[c];
				}
				for(int i=0;i<=b;++i) {
					static mi t[105]; clr(t);
					//j=son<B i=other
					for(int j=0;i+j<=b;++j) for(int x=0;i+j+x<=b;++x) t[j+x]+=w[r][i][j]*pL[j][x];
					clr(w[r][i]);
					for(int j=0;i+j<=b;++j) w[r][i][i+j]=pR[i][j];
					//j=total size
					mul(w[r][i],t,b);
				}
				//i=link to suf
				for(int i=0;i<=b;++i) for(int j=0;j<i;++j) for(int x=i;x<=b;++x) w[r][j][x]+=w[r][i][x]*C[i][j];
			}
			if(k==1) {
				mi z=lf::z[o]*q;
				for(int i=1;i<=n;++i) for(int j=0;j<=i&&j<=bn;++j) f[to[1][1][j]][i]+=w[1][j][i]*z;
				for(int i=1;i<=n;++i) ans[i]+=f[1][i],f[1][i]=0;
				continue;
			}
			for(int t=1;t<=min(k-1,bn);++t) for(int u=1;u<=m;++u) if(ct[u][t]&&mx[u]<k) {
				for(int i=0,x=ct[u][t];i<x;++i) {
					int v=to[u][t][i]; mi z=mi((x-i)&1?MOD-1:1)*C[x][i];
					for(int j=sz[u];j<=n-sz[u]*4;++j) f[v][j]+=f[u][j]*z;
				}
			}
			for(int u=1;u<=m;++u) if(mx[u]<k) {
				static mi c[105][105];
				int up=(n-sz[u])/k; mi z=1;
				for(int i=1;i<=up;++i) {
					z*=h[u][k];
					for(int x=i;x<=up;++x) for(int y=0,ry=min(x,(up-x)/4);y<=ry;++y) c[x][y]+=z*w[i][y][x];
				}
				for(int x=1;x<=up;++x) for(int y=0,ry=min(x,(up-x)/4);y<=ry;++y) {
					for(int i=sz[u],ri=n-max(sz[u]*4,x*k+y*4*k);i<=ri;++i) {
						g[y][u][i+x*k]+=f[u][i]*c[x][y];
					}
					c[x][y]=0;
				}
				for(int i=sz[u];i<=n-sz[u]*4;++i) g[0][u][i]+=f[u][i],f[u][i]=0;
			}
			for(int c=0;c<=n/k/5;++c) {
				for(int t=min(k-1,bn);t;--t) for(int u=1;u<=m;++u) if(ct[u][t]&&mx[u]<k){
					int e=ct[u][t],L=sz[u]+c*k,R=c*k*4;
					for(int i=t+1;i<=min(k-1,bn);++i) R+=ct[u][i]*i*4;
					for(int x=0;x<=e;++x) {
						int up=n-R-x*t*4,v=to[u][t][x];
						if(x==e) for(int i=max(up+1,L);i<=n;++i) g[c][u][i]=0;
						else for(int i=L;i<=up;++i) g[c][v][i]+=g[c][u][i]*C[e][x];
					}
				}
				for(int u=1;u<=m&&sz[u]+c*k<=bn;++u) if(mx[u]<k) {
					int v=c?to[u][k][c]:u;
					for(int i=sz[v];i<=n-sz[v]*4;++i) f[v][i]+=g[c][u][i],g[c][u][i]=0;
				}
			}
			for(int i=1;i<=n;++i) ans[i]+=f[1][i],f[1][i]=0;
		}
	}
	ans[1]=ans[2]=ans[3]=0;
	for(int i:{0,1,2,3}) ans[lf::sz[i]]+=pq[__builtin_popcount(lf::tr[i])]*lf::w[i];
	for(int i=1;i<=n;++i) cout<<(ans[i]*fac[i-1]).x<<"\n";
	return 0;
}