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

题目

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

思路

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