Wannafly挑战赛23F-计数【原根,矩阵树定理,拉格朗日插值】

1|0正题

题目链接:https://ac.nowcoder.com/acm/contest/161/F

---

1|1题目大意

给出$n$个点的一张图，求它的所有生成树中权值和为$k$的倍数的个数。输出答案对$p$取模
$1\leq n,k\leq 100,1\leq m\leq 10^4,p\in[2,10^9]\cup Pri$
数据保证$k\equiv 1(mod\ p)$

---

1|2解题思路

一个想法是把一条边权看做$x^w$的多项式，用矩阵树定理乘起来后$k$的倍数的系数和就是答案。
但是这样系数是$nk$个，显然搞不定。

类似于CF917D-StrangerTree的做法我们可以带入若干个值然后跑矩阵数之后求出若干个点值。

但是这里是$k$的倍数，我们要模拟卷积，可以带入$k$个$g$满足$g^k=1$的就可以了。

这里保证了$k\equiv 1(mod\ p)$，所以我们求出$p$的原根$g$然后带入$g^{\frac{p-1}{k}\times i}(i\in[0,k-1])$就可以了。

然后直接拉插求出$x=0$的点值就好了。

时间复杂度$O(n^3k)$

---

1|3code

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<vector>
#define ll long long
using namespace std;
const ll N=110;
struct node{
	ll x,y,w;
}e[N*N];
ll n,m,k,P,g,x[N],y[N];
vector<ll> p;
ll power(ll x,ll b){
	ll ans=1;
	while(b){
		if(b&1)ans=ans*x%P;
		x=x*x%P;b>>=1;
	}
	return ans;
}
void Prime(){
	ll x=P-1;
	for(ll i=2;i*i<=x;i++)
		if(x%i==0){
			p.push_back(i);
			while(x%i==0)x/=i;
		}
	if(x>1)p.push_back(x);
}
bool check(ll x){
	for(ll i=0;i<p.size();i++)
		if(power(x,(P-1)/p[i])==1)return 0;
	return 1;
}
namespace Matrix{
	ll a[N][N];
	ll det(ll v){
		memset(a,0,sizeof(a));
		ll ans=1;
		for(ll i=1;i<=m;i++){
			ll x=e[i].x,y=e[i].y,w=power(v,e[i].w); 
			(a[x][y]+=P-w)%=P;(a[y][x]+=P-w)%=P;
			(a[x][x]+=w)%=P;(a[y][y]+=w)%=P;
		}
		ll f=0;
		for(ll i=1;i<n;i++){
			for(ll j=i;j<n;j++)
				if(a[j][i]){
					if(i==j)break;
					swap(a[i],a[j]);
					f^=1;break;
				}
			ans=ans*a[i][i]%P;
			ll inv=power(a[i][i],P-2);
			for(ll j=i;j<n;j++)a[i][j]=a[i][j]*inv%P;
			for(ll j=i+1;j<n;j++){
				ll rate=P-a[j][i];
				for(ll k=i;k<n;k++)
					(a[j][k]+=rate*a[i][k])%=P;
			}
		}
		return f?(P-ans):ans;
	}
}
signed main()
{
	scanf("%lld%lld%lld%lld",&n,&m,&k,&P);
	Prime();g=1;
	while(!check(g))
		g++;
	for(ll i=1;i<=m;i++)
		scanf("%lld%lld%lld",&e[i].x,&e[i].y,&e[i].w);
	for(ll i=1;i<=k;i++){
		x[i]=power(g,(P-1)/k*(i-1));
		y[i]=Matrix::det(x[i]);
	}
	ll ans=0;
	for(ll i=1;i<=k;i++){
		ll tmp=1;
		for(ll j=1;j<=k;j++)
			if(i!=j)tmp=tmp*(P-x[j])%P*power(x[i]-x[j],P-2)%P;
		(ans+=tmp*y[i]%P)%=P;
	}
	printf("%lld\n",(ans+P)%P);
	return 0;
}

__EOF__

本文作者：QuantAsk
本文链接：https://www.cnblogs.com/QuantAsk/p/14449695.html
关于博主：退役OIer，GD划水选手
版权声明：本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！
声援博主：如果您觉得文章对您有帮助，可以点击文章右下角【推荐】一下。您的鼓励是博主的最大动力！