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

题面

数数

题解

tm怎么天天数数（

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

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

考虑矩阵树定理。

先推出基尔霍夫矩阵 \(K\)：

\[\begin{aligned} K=&D-A\\ =& \begin{vmatrix} (n-1)x&0&\cdots&0\\ 0&(n-2)x+1&\cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&n-1 \end{vmatrix}_{n\times n} -\begin{vmatrix} 0&x&\cdots&x\\ 1&0&\cdots&x\\ \vdots&\vdots&\ddots&\vdots\\ 1&1&\cdots&0 \end{vmatrix}_{n\times n}\\ =& \begin{vmatrix} (n-1)x&-x&\cdots&-x\\ -1&(n-2)x+1&\cdots&-x\\ \vdots&\vdots&\ddots&\vdots\\ -1&-1&\cdots&n-1 \end{vmatrix}_{n\times n} \end{aligned} \]

其中：

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

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

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

\[\begin{aligned} &[x^m] \begin{vmatrix} (n-1)x&-x&\cdots&-x\\ -1&(n-2)x+1&\cdots&-x\\ \vdots&\vdots&\ddots&\vdots\\ -1&-1&\cdots&n-1 \end{vmatrix}_{(n-1)\times (n-1)}\\ =&(-1)^{n+1}[x^m] \begin{vmatrix} 1&1&\cdots&1&1\\ (n-1)x&-x&\cdots&-x&0\\ -1&(n-2)x+1&\cdots&-x&0\\ \vdots&\vdots&\ddots&\vdots&\vdots\\ -1&-1&\cdots&n-1&0 \end{vmatrix}_{n\times n}\\ \end{aligned} \]

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

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

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

\[\begin{aligned} =&(-1)^{n+1}[x^m] \begin{vmatrix} 1&1&\cdots&1&1\\ (n-1)x+1&1-x&\cdots&1-x&1\\ 0&(n-2)x+2&\cdots&1-x&1\\ \vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&\cdots&n&1 \end{vmatrix}_{n\times n}\\ \end{aligned} \]

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

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

设 \(C_i\) 表示由 \(T\) 矩阵第 \(i\sim n\) 列和第 \(i\sim n\) 行构成的矩阵，\(\det\nolimits_i\) 表示 \(C_i\) 的行列式。

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

\[\begin{aligned} \det\nolimits_i=&(1-x)\det\nolimits_{i+1}-\big[(n-a_i)x+a_i\big]\det\nolimits_{i+1}\\ =&\big[(a_i-n-1)x+1-a_i\big]\det\nolimits_{i+1} \end{aligned} \]

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

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

但你发现 \(\det\nolimits_1\) 好像不太好处理。

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

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

    在这里插入图片描述

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

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

    在这里插入图片描述

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

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

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

    在这里插入图片描述

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

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

    在这里插入图片描述

    那么对于这种情况，显然它的贡献为：

    \[\begin{aligned} &(-1)^{i-1}\left(\prod_{j=1}^{i-1}(n-a_j)x+a_j\right)\det\nolimits_{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{aligned} \]

    设 \(p_j=(a_j-n)x-a_j\)，\(p_{l..r}=\prod\limits_{j=l}^r p_j\)，\(f_{l..r}=\prod\limits_{j=l}^r f_j\)。

    然后枚举每一个 \(i\)，统计第一列选第二个时的答案：

    \[\sum_{i=2}^n p_{1..i-1}f_{i+1..n} \]

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

    具体来说，假设当前要求 \(S_{l,r}=\sum\limits_{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\nolimits_{1}=\det\nolimits_2+S_{2,n} \]

输出 \(\det_1\) 的 \(m\) 次项系数即可。

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

注意所有的东西都可以用一次分治 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
*/
posted @ 2022-10-30 14:02  ez_lcw  阅读(55)  评论(0)  编辑  收藏  举报