【牛客Wannafly挑战赛23 F】计数

题目

题目链接：https://ac.nowcoder.com/acm/contest/161/F
小O有一个n个点，m条边的边带权无向图。小O希望从这m条边中，选出一些边，使得这些边能构成这n个点的生成树。但他还有个幸运数字k。因此他希望最终选出来的这些边的权值和是k的倍数。他想知道最终有多少种可能的方案选出合法的生成树。答案可能很大，幸好小O还有一个幸运质数p。你只需要输出答案对p取模即可。
\(n,k\leq 100,m\leq 10000,p\leq 10^9,p\bmod k=1\)。

思路

生成树计数问题再加上点数很小，考虑矩阵树定理。但是我们要求的是边权和为 \(k\) 倍数不好搞，套路性把每条边边权设为 \(x^d\)，其中 \(d\) 是这条路的长度，\(x\) 是常数。那么我们只需要知道跑完矩阵树定理之后多项式 \(x^{ik}(i\in \text{Z})\) 的系数之和即可。
但是多项式项数太多无法直接求解，由于需要求的是 \(k\) 的倍数，考虑选用恰当的 \(x\) 使得多项式对 \(x^k\) 取模。这样常数项的系数就是答案。
那么显然我们选取的 \(x\) 需要满足 \(x^k\equiv 1\pmod p\)。注意到 \(p\) 是质数且 \(p\bmod k=1\)，等价于 \(k|\varphi(p)\)；设 \(g\) 是 \(p\) 的原根，又因为 \(g^{\varphi(p)}\equiv 1\pmod p\)，我们可以取 \(x=g^{i\cdot \frac{p-1}{k}}(i\in \text{Z}∩[0,k)\ )\)，这样 \(x^k=g^{i\varphi(p)}\equiv 1\pmod p\) 了。
所以我们可以分别带入 \(k\) 个值，可以得到这个 \(k-1\) 次多项式的 \(k\) 个点值，求常数项就直接拉格朗日插值带入 \(0\) 即可。
时间复杂度 \(O(kn^3+km\log p)\)。

代码

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

const int N=110,M=10010;
int n,m,t,G,MOD,U[M],V[M],D[M];
ll ans,f[N][N],x[N],y[N];

ll fpow(ll x,ll k)
{
	ll ans=1;
	for (;k;k>>=1,x=x*x%MOD)
		if (k&1) ans=ans*x%MOD;
	return ans;
}

int findg(int p)
{
	vector<int> d;
	int k=p-1;
	for (int i=2;i*i<=k;i++)
		if (k%i==0)
		{
			d.push_back(i);
			while (k%i==0) k/=i;
		}
	if (k>1) d.push_back(k);
	for (int i=2;;i++)
	{
		bool flag=1;
		for (int j=0;j<d.size();j++)
			if (fpow(i,(p-1)/d[j])==1) { flag=0; break; }
		if (flag) return i;
	}
	return 123456;
}

ll det()
{
	ll res=1;
	for (int i=1;i<n;i++)
	{
		for (int j=i;j<n;j++)
			if (f[i][j])
			{
				if (i!=j) res=-res;
				for (int k=1;k<n;k++)
					swap(f[i][k],f[j][k]);
				break;
			}
		res=res*f[i][i]%MOD;
		ll inv=fpow(f[i][i],MOD-2);
		for (int j=i+1;j<n;j++)
			if (f[j][i])
			{
				ll base=f[j][i]*inv%MOD;
				for (int k=1;k<n;k++)
					f[j][k]=(f[j][k]-f[i][k]*base)%MOD;
			}
	}
	return res;
}

int main()
{
	scanf("%d%d%d%d",&n,&m,&t,&MOD);
	for (int i=1;i<=m;i++)
		scanf("%d%d%d",&U[i],&V[i],&D[i]);
	G=findg(MOD);
	for (int j=0;j<t;j++)
	{
		memset(f,0,sizeof(f));
		x[j]=fpow(G,(MOD-1)/t*j);
		for (int i=1;i<=m;i++)
		{
			ll pw=fpow(x[j],D[i]);
			f[U[i]][V[i]]=(f[U[i]][V[i]]-pw)%MOD;
			f[V[i]][U[i]]=(f[V[i]][U[i]]-pw)%MOD;
			f[U[i]][U[i]]=(f[U[i]][U[i]]+pw)%MOD;
			f[V[i]][V[i]]=(f[V[i]][V[i]]+pw)%MOD;
		}
		y[j]=det();
	}
	for (int i=0;i<t;i++)
	{
		ll mul=1;
		for (int j=0;j<t;j++)
			if (i!=j) mul=mul*(-x[j])%MOD*fpow(x[i]-x[j],MOD-2)%MOD;
		ans=(ans+mul*y[i])%MOD;
	}
	printf("%lld",(ans%MOD+MOD)%MOD);
	return 0;
}