【XSY3918】数数（矩阵树定理，分治NTT）

题面

数数

题解

tm怎么天天数数（

考虑一个点数为 $n$ 的有向完全图：对于边 $i\to j$ 来说，若 $i<j$ 则边权为 $x$，否则边权为 $1$。

显然答案就是这张图中以 $r$ 为根的且边权乘积为 $x^m$ 的内向树个数。

考虑矩阵树定理。

先推出基尔霍夫矩阵 $K$：

$$\begin{array}{rl}K=& D-A\\ =& {\left|\begin{array}{cccc}(n-1)x& 0& \cdots & 0\\ 0& (n-2)x+1& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & n-1\end{array}\right|}_{n\times n}-{\left|\begin{array}{cccc}0& x& \cdots & x\\ 1& 0& \cdots & x\\ ⋮& ⋮& \ddots & ⋮\\ 1& 1& \cdots & 0\end{array}\right|}_{n\times n}\\ =& {\left|\begin{array}{cccc}(n-1)x& -x& \cdots & -x\\ -1& (n-2)x+1& \cdots & -x\\ ⋮& ⋮& \ddots & ⋮\\ -1& -1& \cdots & n-1\end{array}\right|}_{n\times n}\end{array}$$

其中：

• $D$ 为度数矩阵（加权），只有 $(i,i)$ 有值。同时这是内向树，统计的是每个点的加权出度。

• $A$ 为邻接矩阵（加权），除了 $(i,i)$ 都有值。

由于需要以 $r$ 为根，所以我们去掉第 $r$ 行和第 $r$ 列，变成一个 $(n-1)\times (n-1)$ 的矩阵（设去掉第 $r$ 行和去掉第 $r$ 列之前的 $n\times n$ 的矩阵为 $S$）。那么答案为：

$$\begin{array}{rl}& [x^m]{\left|\begin{array}{cccc}(n-1)x& -x& \cdots & -x\\ -1& (n-2)x+1& \cdots & -x\\ ⋮& ⋮& \ddots & ⋮\\ -1& -1& \cdots & n-1\end{array}\right|}_{(n-1)\times (n-1)}\\ =& (-1{)}^{n+1}[x^m]{\left|\begin{array}{ccccc}1& 1& \cdots & 1& 1\\ (n-1)x& -x& \cdots & -x& 0\\ -1& (n-2)x+1& \cdots & -x& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ -1& -1& \cdots & n-1& 0\end{array}\right|}_{n\times n}\end{array}$$

上面这里我们是新加了一行和一列，我们来证明这个等式为什么是对的：

不妨设原来的行列式的值为 $v$，新加了一行一列的行列式的值为 $v^{\prime }$。结合行列式的定义，显然当矩阵第一行选最后一个 $1$ 的时候才对行列式的值有贡献，此时会比没有新添加一行一列时多出 $n-1$ 个逆序对。所以得到 $v\times (-1{)}^{n-1}=v^{\prime }$，也就得到 $v=v^{\prime }\times (-1{)}^{n-1}$。

那么由于矩阵的一行加上另一行的倍数后行列式不变，所以不妨将第 $2\sim n$ 行都加上第 $1$ 行，那么此时原式变为：

$$\begin{array}{rl}=& (-1{)}^{n+1}[x^m]{\left|\begin{array}{ccccc}1& 1& \cdots & 1& 1\\ (n-1)x+1& 1-x& \cdots & 1-x& 1\\ 0& (n-2)x+2& \cdots & 1-x& 1\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & n& 1\end{array}\right|}_{n\times n}\end{array}$$

设最后推出来的这个矩阵为 $T$。（注意：$T$ 不是上三角矩阵，因为对角线是 $1,\underset{(n-2)\text{个}}{\underbrace{1-x,\cdots ,1-x}},1$）

由于我们一开始删掉了第 $r$ 行和第 $r$ 列，所以我们不妨设 $T$ 矩阵的第 $i$ 列代表的是一开始 $S$ 矩阵的第 $a_i$ 列，那么易知 $a_i=i+[i\ge r](1\le i<n)$。

设 $C_i$ 表示由 $T$ 矩阵第 $i\sim n$ 列和第 $i\sim n$ 行构成的矩阵，${\det }_i$ 表示 $C_i$ 的行列式。

按第一列选哪个分类讨论，易得下面的递推式：（$2\le i<n$，${\det }_n=1$）

$$\begin{array}{rl}{\det }_i=& (1-x){\det }_{i+1}-[(n-a_i)x+a_i]{\det }_{i+1}\\ =& [(a_i-n-1)x+1-a_i]{\det }_{i+1}\end{array}$$

设 $f_i=(a_i-n-1)x+1-a_i$（特判 $f_n=1$），显然有 ${\det }_2=\prod _{i=2}^n{f}_i$。

这是简单的分治 NTT，可以轻松实现。

但你发现 ${\det }_1$ 好像不太好处理。

我们还是按第一列哪一个分类讨论：

• 如图，如果我们第一列选了第一个数 $1$，那么我们之后选的数只能在红框里面选了：

  

  又由于第一列第一个数不会和后面的数产生逆序对，所以易知这种情况的贡献为 $1\times f_2=f_2$。

• 如图，如果我们第一列选了第二个数 $(n-a_1)x+a_1$：

  

  显然，我们并不能用 $f_2$ 转移，因为第一行后面是 $1$，而第二行后面是 $1-x$，它们不相同。

  那么我们枚举第一行选了哪一个，设第一行选了第 $i$ 列的 $1$：（$1<i\le n$）

  （十分抱歉，笔者下面的所有图中所有的 $(n-a_i)x+a_i$ 都写成了 $(n-a_i)+a_i$，请自行脑补一个 $x$ 上去）

  

  然后你发现，第二列只能选第三个数，因为前两行都被选过了，而选第四行以后都是 $0$ 没有贡献。

  然后再推一推就能发现，对于前 $i-1$ 列中的每一列 $j(1\le j<i)$，都只能选第 $j+1$ 行，而对于第 $i$ 列后面的所有列，它们可以且只能在 $C_{i+1}$ 里面选：

  

  那么对于这种情况，显然它的贡献为：
  $$\begin{array}{rl}& (-1{)}^{i-1}\left(\prod _{j=1}^{i-1}(n-a_j)x+a_j\right){\det }_{i+1}\\ =& \left(\prod _{j=1}^{i-1}(a_j-n)x-a_j\right)\left(\prod _{j=i+1}^n{f}_i\right)\end{array}$$
  设 $p_j=(a_j-n)x-a_j$，$p_{l..r}=\prod _{j=l}^r{p}_j$，$f_{l..r}=\prod _{j=l}^r{f}_j$。

  然后枚举每一个 $i$，统计第一列选第二个时的答案：
  $$\sum _{i=2}^n{p}_{1..i-1}{f}_{i+1..n}$$
  这是简单的分治 NTT，可以轻松实现。

  具体来说，假设当前要求 $S_{l,r}=\sum _{i=l}^r{p}_{l-1..i-1}{f}_{i+1...r+1}$，而且左右都分治好了，那么：
  $$S_{l,r}=S_{l,mid}{f}_{mid+2..r+1}+p_{l-1..mid-1}{S}_{mid+1,r}$$

那么两种情况都处理完了，最后得到：
$${\det }_1={\det }_2+S_{2,n}$$
输出 ${\det }_1$ 的 $m$ 次项系数即可。

时间复杂度 $O(n{\log }^{2}n)$。

注意所有的东西都可以用一次分治 NTT 求出来，不需要求两次。

代码如下：

#include<bits/stdc++.h>

#define LN 18
#define N 100010
#define mod 998244353

using namespace std;

inline int read()
{
	int x=0,f=1;
	char ch=getchar();
	while(ch<'0'||ch>'9')
	{
		if(ch=='-') f=-1;
		ch=getchar();
	}
	while(ch>='0'&&ch<='9')
	{
		x=(x<<1)+(x<<3)+(ch^'0');
		ch=getchar();
	}
	return x*f;
}

namespace modular
{
	int add(int x,int y){if((x+=y)>=mod)x-=mod; return x;}
	int dec(int x,int y){if((x-=y)<0)x+=mod; return x;}
	int mul(int x,int y){return 1ll*x*y%mod;}
}

using namespace modular;

int n,rt,m;
int rev[N<<1],w[LN][N<<1][2];
int sl[N<<1],sr[N<<1],fl[N<<1],fr[N<<1],pl[N<<1],pr[N<<1];
int ns[N<<1],nf[N<<1],np[N<<1];

vector<int>f[N<<2],p[N<<2],s[N<<2];

int poww(int a,int b)
{
	int ans=1;
	while(b)
	{
		if(b&1) ans=mul(ans,a);
		a=mul(a,a);
		b>>=1;
	}
	return ans;
}

void init(int limit)
{
	for(int bit=0,mid=1;mid<limit;bit++,mid<<=1)
	{
		int len=mid<<1;
		int gn=poww(3,(mod-1)/len);
		int ign=poww(gn,mod-2);
		int g=1,ig=1;
		for(int j=0;j<mid;g=mul(g,gn),ig=mul(ig,ign),j++)
			w[bit][j][0]=g,w[bit][j][1]=ig;
	}
}

void NTT(int *a,int limit,int opt)
{
	opt=(opt<0);
	for(int i=0;i<limit;i++)
		rev[i]=(rev[i>>1]>>1)|((i&1)*(limit>>1));
	for(int i=0;i<limit;i++)
		if(i<rev[i]) swap(a[i],a[rev[i]]);
	for(int bit=0,mid=1;mid<limit;bit++,mid<<=1)
	{
		for(int len=mid<<1,i=0;i<limit;i+=len)
		{
			for(int j=0;j<mid;j++)
			{
				int x=a[i+j],y=mul(w[bit][j][opt],a[i+mid+j]);
				a[i+j]=add(x,y),a[i+mid+j]=dec(x,y);
			}
		}
	}
	if(opt)
	{
		int tmp=poww(limit,mod-2);
		for(int i=0;i<limit;i++)
			a[i]=mul(a[i],tmp);
	}
}

void solve(int k,int l,int r)
{
	if(l==r)
	{
		if(l<n-1)
		{
			int a=(l+1)+((l+1)>=rt);
			f[k].push_back(dec(1,a));
			f[k].push_back(dec(a,n+1));
		}
		else
		{
			f[k].push_back(1);
			f[k].push_back(0);
		}
		int a=(l-1)+((l-1)>=rt);
		p[k].push_back(dec(0,a));
		p[k].push_back(dec(a,n));
		s[k].push_back(mul(f[k][0],p[k][0]));
		s[k].push_back(add(mul(f[k][1],p[k][0]),mul(f[k][0],p[k][1])));
		s[k].push_back(mul(f[k][1],p[k][1]));
		return;
	}
	int mid=(l+r)>>1,lc=k<<1,rc=k<<1|1;
	solve(lc,l,mid),solve(rc,mid+1,r);
	
	for(int i=0,size=s[lc].size();i<size;i++) sl[i]=s[lc][i];
	for(int i=0,size=s[rc].size();i<size;i++) sr[i]=s[rc][i];
	for(int i=0,size=f[lc].size();i<size;i++) fl[i]=f[lc][i];
	for(int i=0,size=f[rc].size();i<size;i++) fr[i]=f[rc][i];
	for(int i=0,size=p[lc].size();i<size;i++) pl[i]=p[lc][i];
	for(int i=0,size=p[rc].size();i<size;i++) pr[i]=p[rc][i];
	
	int limit=1;
	while(limit<=(r-l+2)) limit<<=1;
	NTT(sl,limit,1),NTT(sr,limit,1),NTT(fl,limit,1),NTT(fr,limit,1),NTT(pl,limit,1),NTT(pr,limit,1);
	for(int i=0;i<limit;i++)
	{
		ns[i]=add(mul(sl[i],fr[i]),mul(pl[i],sr[i]));
		nf[i]=mul(fl[i],fr[i]);
		np[i]=mul(pl[i],pr[i]);
	}
	NTT(ns,limit,-1),NTT(nf,limit,-1),NTT(np,limit,-1);
	
	for(int i=0;i<=(r-l+2);i++) s[k].push_back(ns[i]);
	for(int i=0;i<=(r-l+1);i++) f[k].push_back(nf[i]);
	for(int i=0;i<=(r-l+1);i++) p[k].push_back(np[i]);
	
	for(int i=0;i<limit;i++) sl[i]=sr[i]=fl[i]=fr[i]=pl[i]=pr[i]=ns[i]=nf[i]=np[i]=0;
}

int main()
{
	n=read(),rt=read(),m=read();
	if(n==1)
	{
		cout<<1<<endl;
		return 0;
	}
	int limit=1;
	while(limit<=n) limit<<=1;
	init(limit);
	solve(1,2,n);
	if(2<n-1)
	{
		int a=2+(2>=rt);
		fl[0]=dec(1,a),fl[1]=dec(a,n+1);
	}
	else fl[0]=1;
	for(int i=0,size=f[1].size();i<size;i++) fr[i]=f[1][i];
	NTT(fl,limit,1),NTT(fr,limit,1);
	for(int i=0;i<limit;i++)
		fl[i]=mul(fl[i],fr[i]);
	NTT(fl,limit,-1);
	int ans=add(m<n-1?fl[m]:0,s[1][m]);
	printf("%d\n",(n&1)?ans:dec(0,ans));
	return 0;
}
/*
4 3 2
*/