这题超级困难吧，感觉不止紫啊。

简化题意：

给你一颗树，保证每个点的父亲编号小于自身。

定义一次任务为给定集合点 $x$，然后选定 $S$ 使得 $\text{LCA}(S)=x$.

要求每个点至多出现在一个 $S$ 中，求有多少种完成 $k$ 个任务的不同方案。

$n\le 10^6,k\le 2000$，保证 $x_i$ 两两不同。

---

在讲这个题的做法之前，先讲一个思考方向类似的题。

前置：P13525 [KOI 2025 #2] 新的情缘。

一道很神奇的题，感觉肯定有蓝。

首先想一想题目保证区间无交有什么用，考虑若干个不交区间组的后缀，这个后缀的男孩子和女孩子的数量相等。如果前面一个男孩子选了这个后缀中的女孩子，则这个后缀中的某个男孩子因为选不到前面的女孩子而失配。

因此把区间包含关系建成一个森林，则每颗小树互相独立，总方案数为每颗小树的方案之积。

先不管不能复合的限制，统计整颗树的方案数。认为女孩子有择偶权，考虑括号匹配（保证每个女孩子都有男孩子可选），则每个右括号加入栈的时刻，这个女孩子可以在栈里面选一个男孩子，方案数为栈的大小（从带走最后一个左括号变成挑一个左括号带走）。进一步地，这个右括号入栈时栈里面的左括号个数，其实就是包含该区间的区间数量 $+1$，也就是这个右括号对应结点的深度（根的深度为 $1$）。

也就是说，若不考虑不能复合的限制，则 $\text{Ans}_T=\prod \text{dep}_u$.

加上不能复合的限制，考虑容斥，若我们钦定 $|S|$ 中的情侣复合，则这部分的贡献为 $(-1)^{|S|}\prod_{u\not\in S} (\text{dep}_u-\text{cnt}_u)$，其中 $\text{cnt}_u$ 表示 $u$ 祖先链上的钦定复合的情侣数，由括号匹配的过程不难理解。容易想到，令：

$$w_S(u)=\begin{cases} \text{dep}_u-\text{cnt}_u, u\not\in S\\ -1,u\in S \end{cases}$$

于是最终答案可以表示为 $\text{Ans}=\sum_{S\sub [1,\dots n]}\prod_u w_S(u)$.

考虑 dp 算这个玩意，设 $f_{i,j}$ 表示我们钦定 $\text{fa}_i$ 到根链上有 $j$ 对情侣 没有钦定复合，$i$ 子树的上面那个式子的值。转移是简单的，复杂度 $\sum_u \text{deg}_u\times \text{dep}_u=\sum_{v}\text{dep}_v=O(n^2)$.

End.

---

回到本题，每个点只能在一个集合中出现的条件非常强，我们不得不同时考虑所有任务。定义给定的集合点为 任务点，则每个点都可以选择完成它到根链上的一个任务点，或者摆烂。此时如果不要求任务完成，一次 dfs 就可以算得方案数。一个任务被完成，当且仅当它所在结点选择完成它，或者至少两个儿子的子树内的点选择完成它。

任务都被完成也太难刻画了，考虑容斥，仿照前置题那样，设 $f_{i,j}$ 表示 $i$ 子树，钦定 $\text{fa}_i$ 到根链上有 $j$ 个任务 没有被钦定不满足。

若 $u$ 上没有任务，则转移就是 $f_{u,j}=(j+1)\prod_{v\in \text{son}_u}f_{v,j}$.

若 $u$ 上有任务，转移比较复杂：

• 若 $u$ 的任务由 $u$ 直接完成，$f_{u,j}\gets 1\cdot \prod_k f_{k,j+1}$.
• 若 $u$ 的任务没有由 $u$ 完成，辅助转移出 $g_{0/1/2}$ 表示 $i$ 子树钦定值为 $j$，选择了 $0/1/2$ 个有效儿子（显然 $\ge 2$ 与 $2$ 等价）的加权贡献和，辅助转移可以依靠背包完成。
  • 于是，若我们钦定 $u$ 的任务不满足，$f_{i,j}\gets (j+1)\cdot(-g_{0/1})$.
  • 否则，任务是否满足是随意的，$f_{i,j}\gets (j+1)\cdot g_{0/1/2}=(j+1)\prod_{v}f_{v,j+1}$.
  • 这两者加和其实也就是把 $(j+1)\cdot g_{2}$ 加入答案。可以认为 这一篇题解 就是把 $g_0$ 和 $g_1$ 直接展开了。

其中，“有效儿子”被认为是 $u$ 计入其任务列表且 $u$ 确实被完成的那些儿子，对应加入背包的系数为 $f_{v,j+1}-f_{v,j}$，否则对应那些 $u$ 根本没计入任务列表的儿子，背包系数为 $f_{v,j}$.

$g$ 的转移是：

$$g'_0\gets g_0\times f_{v,j}\\ g'_1\gets g_1\times f_{v,j}+g_0\times (f_{v,j+1}-f_{v,j})\\ g'_2\gets g_2\times f_{v,j+1}+g_1\times (f_{v,j+1}-f_{v,j})$$

复杂度为 $\sum k\cdot \text{deg}_u=O(nk)$.

由于上面写的比较抽象，给代码：$O(nk)$ Code.

因为非任务点远多于任务点，所以大多数转移都是 $u$ 上无任务的简单转移，这里用到了虚状态合并的思想，假设 $u,\text{fa}_u$ 都无任务，则 $p=\text{fa}_u$ 的式子是：

$$\begin{aligned} f_{p,j}&=(j+1)f_{u,j}\prod_{x\in \text{son}_p,x\not=u} f_{x,j}\\ &=(j+1)^2\prod_{y\in son_u}f_{y,j}\prod_{x\in \text{son}_p,x\not=u} f_{x,j} \end{aligned}$$

这个转移和普通无任务的转移根本没有本质区别，只是系数不同。若我们把 $p$ 和 $u$ 两个点合并，然后把 $u$ 的儿子都拉到 $p$ 上，转移基本上是等价的，只需记一个 $c$ 表示这个点是 $c$ 个点叠合而成的就行了。

因此得到一个做法：按深度从大到小考虑每个点 $u$，若 $u$ 和 $\text{fa}_u$ 都是非任务点，则把 $u$ 的所有儿子并入 $\text{fa}_u$ 的邻接表，并在 $\text{fa}_u$ 的邻接表中删除点 $u$，处理完成后，跑前面的 $O(nk)$ 的 dp 算法。

如果你使用 std::set 维护邻接表，加以启发式合并，可以把这部分做到 $O(n\log^2 n)$，但是似乎根本没人卡这个做法，所以我用 std::vector 加上暴力合并 直接过了，非常逆天。

你可能会担心这个做法改变了树结构，毕竟其他做法都是基于虚树的，且保留了任务点的父亲为“骨架点”，以保证关键点不会挪动到某棵错误的子树中。但是，观察转移形式可以发现 非任务点的儿子数量并不重要，即使一个任务点错误地成为了一个非任务点的儿子，直接作乘法得到的结果也是正确的。

这个做法太笨蛋了，考虑缩树的过程本质上就是缩掉所有非关键点之间的连边，利用 并查集 把非关键点合成联通块，然后再连树边，就可以得到 $O(n\alpha(n))$ 建新树。

合并之后，不会存在非任务点之间的连边，因而树会变成任务点形成的虚树并上任务点下属的若干非任务叶子（每个这样的叶子代表着原树上的某一颗极大不含任务点的子树）。但是，叶子数量可能是 $O(n)$ 的，所以直接这样复杂度仍然不对。

考虑到这些叶子的“权值”$c$ 之和恰为 $n-k=O(n)$，因而不同的权值种类数只有 $O(\sqrt n)$，如果我们可以打包处理一种“权值”的所有叶子，总复杂度就可以接受了。

首先，“权值”相同的叶子 $f$ 数组完全相同，有 $f_{i,j}=(j+1)^c$，毕竟我们都把它缩成叶子了嘛。因而小改一下上面的背包过程即可，假设这种叶子有 $w$ 个，转移如下：

$$G_0=(j+1)^c,G_1=(j+2)^c-(j+1)^c,G_2=(j+2)^c\\ g'_0\gets g_0\times G_0^w\\ g'_1\gets g_1\times G_0^w+g_0\times wG_1G_0^{w-1}\\ g'_2\gets g_2\times G_2^w+g_1\times (G_2^w-G_0^w)+g_0\times (G_2^w-G_0^w-wG_1G_0^{w-1})$$

注意向 $g'_2$ 的转移，我们用容斥刻画“至少”。

于是空点转移的复杂度降为 $k\sum_{u\in key}\sqrt{e_u}\le k\sqrt{k\sum e_u}=O(k\sqrt{nk})$，需要 $O(k\sqrt n)$ 预处理光速幂。

总复杂度是 $O(n\alpha(n)+k^2+k\sqrt{nk})$.

为了合并同一个结点上挂着的相同权值叶子部分好写，我的做法依旧笨蛋，复杂度多加上一个 $O(n\log n)$.

:::info[Code]

#include<bits/stdc++.h>
using namespace std;
#define debug(...) fprintf(stderr,__VA_ARGS__)
struct FSI{
	template<typename T>
	FSI& operator >> (T&res){
		res=0;T f=1;char ch=getchar();
		while(!isdigit(ch)){if(ch=='-') f=-1;ch=getchar();}
		while(isdigit(ch)){res=(res*10)+(ch-48);ch=getchar();}
		res*=f;
		return *this;
	}
} scan;
template<int umod>
struct modi{
	static constexpr int mod=umod;
	int val;
	modi(){val=0;}
	modi(int _v){val=_v;}
	friend modi operator + (const modi &a,const modi &b){
		modi res(a.val+b.val);
		res.val=res.val>=mod?res.val-mod:res.val;
		return res;
	}
	friend modi operator - (const modi &a,const modi &b){
		modi res(a.val+mod-b.val);
		res.val=res.val>=mod?res.val-mod:res.val;
		return res;
	}
	friend modi operator - (const modi &a){
		return a.val?modi(mod-a.val):a;
	}
	friend modi operator * (const modi &a,const modi &b){
		return modi(1ll*a.val*b.val%mod);
	}
	void operator += (const modi &a){
		val+=a.val;
		val=val>=mod?val-mod:val;
	}
	void operator *= (const modi &a){
		val=1ll*val*a.val%mod;
	}
	// friend modi qpow(modi a,int b){
	// 	modi res(1);
	// 	while(b){
	// 		if(b&1) res*=a;
	// 		a*=a,b>>=1;
	// 	}
	// 	return res;
	// }
};
const int mod=998244353;
using modint=modi<mod>;
const int N=1e6+5;
int n,k;
bool key[N];
vector<int> G[N];
int fa[N],cnt[N],cc[N];
vector<modint> f[N];
bool is_leaf[N];
namespace quick_power{
	const int B=1000;
	modint pw[2010][B+5],pwB[2010][N/B+5];
	void init(int k,int n){
		for(int i=1;i<=k+2;i++){
			pw[i][0]=1;
			for(int j=1;j<=B;j++) pw[i][j]=pw[i][j-1]*modint(i);
			pwB[i][0]=1;
			for(int j=1,lim=n/B+1;j<=lim;j++) pwB[i][j]=pwB[i][j-1]*pw[i][B];
		}
	}
	modint qpow(int a,int b){
		return pwB[a][b/B]*pw[a][b%B];
	}
}
using quick_power::qpow;
void dfs(int u){
	for(int v:G[u]){
		cnt[v]=cnt[u]+key[v],dfs(v);
	}
	f[u].resize(cnt[u]+1);
	if(!key[u]){
		if(!G[u].empty()){//如果是孤叶子，不作处理，防止复杂度退化.
			for(int i=0;i<=cnt[u];i++){
				f[u][i]=qpow(i+1,cc[u]);
				for(int v:G[u]) f[u][i]*=f[v][i];
			}
		} else is_leaf[u]=true;
	} else {
		vector<int> vec;
		map<int,int> mp;
		typedef pair<int,int> pii;
		vector<pii> leaf;
		for(int v:G[u]){
			if(is_leaf[v]) mp[cc[v]]+=1;
			else vec.emplace_back(v);
		}
		for(auto t:mp) leaf.emplace_back(t.first,t.second);
		for(int i=0;i<cnt[u];i++){
			modint g[3]={1,0,0},mul(1);
			for(int v:vec){
				g[2]=g[2]*f[v][i+1]+g[1]*(f[v][i+1]-f[v][i]);
				g[1]=g[1]*f[v][i]+g[0]*(f[v][i+1]-f[v][i]);
				g[0]=g[0]*f[v][i];
				mul*=f[v][i+1];
			}
			for(auto [c,w]:leaf){
				modint G2=qpow(i+2,c*w),G0=qpow(i+1,c*w);
				modint G1=(qpow(i+2,c)-qpow(i+1,c))*qpow(i+1,c*(w-1))*w;
				g[2]=g[2]*G2+g[1]*(G2-G0)+g[0]*(G2-G0-G1);
				g[1]=g[1]*G0+g[0]*G1;
				g[0]=g[0]*G0;
				mul=mul*G2;
			}
			// assert((g[0]+g[1]+g[2]).val==mul.val);
			f[u][i]=mul+modint(i+1)*g[2];
		}
	}
}
namespace dsu{
	int fa[N],sz[N];
	inline int find(int x){
		while(x!=fa[x]) x=fa[x]=fa[fa[x]];
		return x;
	}
	inline void init(){
		for(int i=1;i<=n;i++) fa[i]=i,sz[i]=1;
	}
	void merge(int u,int f){
		f=find(f),u=find(u);
		if(f!=u) sz[f]+=sz[u],sz[u]=0,fa[u]=f;
	}
}
int main(){
	scan>>n>>k;
	for(int i=1,x;i<=k;i++){
		scan>>x;
		key[x]=true;
	}
	for(int i=2;i<=n;i++) scan>>fa[i];
	dsu::init(),quick_power::init(k,n);
	for(int i=2;i<=n;i++){
		if(!key[i]&&!key[fa[i]]) dsu::merge(i,fa[i]);
	}
	for(int i=2;i<=n;i++){
		if(dsu::fa[i]==i){
			G[dsu::find(fa[i])].emplace_back(i);
			if(!key[i]) cc[i]=dsu::sz[i];
		}
	}
	cnt[1]=key[1],dfs(1);
	printf("%d\n",f[1][0].val);
	return 0;
}

:::