4.9 省选模拟赛 生成树求和 变元矩阵树定理 生成函数 iDFT 插值法

avatar
avatar

有同学在loj上找到了加强版 所以这道题是可以交的。LINK:生成树求和 加强版

对于30分 爆搜 可实际上我爆搜只过了25分 有同学使用按秩合并并茶几的及时剪枝通过了30分。

const int MAXN=45;
int n,m;
struct wy
{
	int x,y,z;
}t[MAXN];
int w[MAXN],f[MAXN];ll ans;
inline int getfather(int x){return x==f[x]?x:getfather(f[x]);}
inline int add(int x,int y)
{
	int cnt=0;
	int p=1;
	while(x||y)
	{
		int w1=x%3;
		int w2=y%3;
		int cc=(w1+w2)%3;
		cnt+=cc*p;
		p=p*3;x/=3;y/=3;
	}
	return cnt;
}
inline void dfs(int x,int v)
{
	if(m-x+1<n-v)return;
	if(v==n)
	{
		int cnt=0;
		rep(1,n,i)f[i]=i;
		rep(1,n-1,i)
		{
			int xx=getfather(t[w[i]].x);
			int yy=getfather(t[w[i]].y);
			if(xx==yy)return;
			f[xx]=yy;
			cnt=add(cnt,t[w[i]].z);
		}
		ans=(ans+cnt)%mod;
		return;
	}
	w[v]=x;
	dfs(x+1,v+1);
	dfs(x+1,v);
}
int main()
{
	freopen("sum.in","r",stdin);
	freopen("sum.out","w",stdout);
	get(n);get(m);
	rep(1,m,i)
	{
		int get(x);int get(y);int get(z);
		t[i]=(wy){x,y,z};
	}
	if(m<=35){dfs(1,1);putl(ans);}
	return 0;
}

之所以这样做是因为虽然这道题很容易想到矩阵树定理,按位处理 但是依然无法做出来。

矩阵树除了能求方案还可以求边上的权值积的和,即变元矩阵树定理 具体理解可以把边权想象成有边权条重边。

但是对于这道题来说还是不行 求的是生成树边权和的和 把积变和只有是在指数上了。

而且还有三进制的加法 非常的不可做 但是可以按位处理。

把边权变到指数上只能是 生成函数了 列出生成函数 每一条边的权值用生成函数表示。

利用高斯消元解行列式 最后求出的多项式就好了。

而直接使用多项式无法方便的求出行列式 中间由于矩阵树定理还无法取模 所以有多项式乘法除法啥的 所以变的更不可做了。

但是考虑最后的次数 其实并不大最多2n次 可以取2n+1个数带入 然后利用插值求解原多项式。

复杂度\(n^4log\)

为了通过加强版 我们考虑一种比较简单的方式 如果能找到\(x^3=1(mod P)\)

那么我么最后的多项式就只有三项且是正确的了 考虑到复平面单位根。

由于最后有三项 所以 把复平面平分三份的复平面单位根即可。

考虑最后的插值 利用复平面单位根的相加相乘操作可以轻松的插值出来。不需要再次高斯消元或者拉格朗日插值。

非常的巧妙。

const ll MAXN=100010,inv2=(mod+1)>>1,inv3=(mod+1)/3,g3=82062379,maxn=110;
inline ll ksm(ll b,ll p,ll cc){ll cnt=1;while(p){if(p&1)cnt=cnt*b%cc;p=p>>1;b=b*b%cc;}return cnt;}
struct wy
{
	ll a,b;
	inline wy(ll aa=0,ll bb=0){a=aa;b=bb;}
	inline wy operator +(wy w){return wy((a+w.a)%mod,(b+w.b)%mod);}
	inline wy operator -(wy w){return wy((a-w.a+mod)%mod,(b-w.b+mod)%mod);}
	inline wy operator *(ll y){return wy(a*y%mod,b*y%mod);}
	inline wy operator *(wy y){return wy((a*y.a%mod-b*y.b%mod+mod)%mod,(a*y.b+b*y.a)%mod);}
	inline wy operator /(ll y)
	{
		ll inv=ksm(y,mod-2,mod);
		return wy(a*inv%mod,b*inv%mod);
	}
	inline wy operator /(wy y)
	{
		return wy(a,b)*wy(y.a,mod-y.b)/((y.a*y.a+y.b*y.b)%mod);
	}
}a[maxn][maxn],cs[3],inv,w[3],invw[3];
int pow3[maxn],eu[maxn*maxn],ev[maxn*maxn],ec[maxn*maxn];
ll n,m;
inline wy det()
{
	wy ret=wy(1,0);
	rep(1,n-1,i)
	{
		ll p=i;
		rep(i,n-1,j)if(a[j][i].a||a[j][i].b){p=j;break;}
		if(i!=p)
		{
			rep(1,n-1,k)swap(a[i][k],a[p][k]);
			ret=ret*(mod-1);
		}
		rep(i+1,n-1,j)
		{
			inv=a[j][i]/a[i][i];
			rep(i,n-1,k)a[j][k]=a[j][k]-inv*a[i][k];
		}
		ret=ret*a[i][i];
	}
	return ret;
}
inline void idft(wy *p)
{
	wy c[3];
	c[0]=p[0];c[1]=p[1];c[2]=p[2];
	p[0]=c[0]+c[1]+c[2];
	p[1]=c[0]+invw[1]*c[1]+c[2]*invw[2];
	p[2]=c[0]+invw[2]*c[1]+c[2]*invw[1];
	rep(0,2,i)p[i]=p[i]*inv3;
}
inline ll calc(ll p)
{
	wy w0,w1,w2,c;
	rep(0,2,i)
	{
		memset(a,0,sizeof(a));
		w0=wy(1,0);
		w1=w[i];
		w2=w[(i+i)%3];
		rep(1,m,j)
		{
			ll k=ec[j]/pow3[p]%3;
			c=(k==1?w1:k==2?w2:w0);
			a[eu[j]][eu[j]]=a[eu[j]][eu[j]]+c;
			a[ev[j]][ev[j]]=a[ev[j]][ev[j]]+c;
			a[eu[j]][ev[j]]=a[eu[j]][ev[j]]-c;
			a[ev[j]][eu[j]]=a[ev[j]][eu[j]]-c;
		}
		cs[i]=det();
	}
	idft(cs);
	return (cs[1].a+cs[2].a*2)*pow3[p]%mod;
}
signed main()
{
	freopen("1.in","r",stdin);
	//freopen("sum.out","w",stdout);
	get(n);get(m);ll ans=0;
	w[0]=wy(1,0);w[1]=wy(mod-inv2,g3*inv2%mod);
	w[2]=wy(mod-inv2,mod-g3*inv2%mod);
	invw[0]=w[0]/w[0];
	invw[1]=w[0]/w[1];
	invw[2]=w[0]/w[2];
	rep(1,m,i)get(eu[i]),get(ev[i]),get(ec[i]);
	pow3[0]=1;
	rep(1,10,i)pow3[i]=pow3[i-1]*3;
	rep(0,10,i)ans=(ans+calc(i))%mod;
	putl((ans+mod)%mod);return 0;
}
posted @ 2020-04-11 16:23  chdy  阅读(158)  评论(0编辑  收藏  举报