洛谷 P5933 [清华集训2012]串珠子 题解

原题：[清华集训2012]串珠子

加强版：串珠子 加强版

这道题基本思路与之前写过题解的轰炸方案一样，可见这是一种常用的套路了。

大概翻了翻题解，发现大部分都是用的状压dp+枚举子集，时间复杂度 \(O(3^n)\)。

这里提供一个 \(O(n^2 2^n)\) 的做法，实测可以在 2s 内通过 \(n<=20\) 的数据，目前排全站第一。

前置知识：快速莫比乌斯变换、子集卷积

前面的推式子都差不多，但这里还是说一下。

第一步，推式子

我们定义：

\(f_S\) 表示\(S\)集合构成连通图的方案数；

\(g_S\) 表示\(S\)集合乱连的总方案（即不保证连通）。

我们通过递推的方式求 \(g\) ：

\[g_S=g_{S'} \prod\limits_{q\in S'} (c_{p,q}+1) \]

其中 \(S'\) 为 \(S\) 中去掉任意某元素 \(p\) 所得，“ \(+1\) ”即可以不连。

直接求连通方案 \(f\) 不好求，但我们考虑算不连通方案，再用总方案去减。

对于不连通方案，只要图有两个以上的连通块，这个图就不连通了。我们抓住这个特点，考虑枚举其中一个连通块，其他的点任其自生自灭，只要不连向我们枚举的连通块，随便连，这样的图肯定是不连通的。

但直接枚举可能算重，即我们现在枚举了一个连通块 \(T\)，之后枚举另一个连通块 \(T'\) 时，\(T\) 出现在了自生自灭的那些点中，这显然会重复计算。

只要我们先后枚举的连通块总是有交，就一定不会算重了。考虑从 \(S\) 中抽出一个点 \(p\)，枚举的所有连通块中均包含 \(p\)，就不会出现上述情况。跟上个式子一样，\(S'\) 为 \(S\) 中去掉任意某元素 \(p\) 所得：

\[f_S=g_S-\sum\limits_{T\subset S'}g_{S'-T}f_{T\cup \{p\}} \]

这里图方便用减号代替了补集符号，大家看的懂就行。

做到这一步，已经可以用 \(O(3^n)\) 枚举子集做了。下面我将对式子做一些变换，来做到更优复杂度。

第二步，一点点变换

看着这个式子，是不是挺像子集卷积？那我们就把它变得更子集卷积。

注意到原式子中 \(T\) 是不能取 \(S'\) 的，于是我们移项并合并：

\[\sum\limits _{T\subseteq S'}g_{S'-T}f_{T\cup \{p\}} =g_S \]

其实式子已经变完了（真·一点点变换），但为了求出答案，我们改变一下 \(p\) 的定义。

\(p\) 可以是 \(S\) 中的任意元素，于是我们钦定 \(p\) 为 \(S\) 中最高位元素，即 \(p=highbit(S)\)。

设总点数为 \(n\)，考虑 \(S\) 取遍 \([2^{n-1},2^n-1]\)，则 \(p=2^{n-1}\)，\(S'\in [0,2^{n-1}-1]\) 。

回头看我们的式子，发现其实就是 \(g\) 前半段与 \(f\) 后半段子集卷积后变成了 \(g\) 后半段。

这里设 \(G_0\) 为 \(g\) 前半段，\(G_1\) 为后半段，\(f\) 同理：

\[G_0*F_1=G_1 \]

\[F_1=\frac{G1}{G0} \]

其中 \(*\) 为子集卷积，除号即为其逆运算。

而 \(g\) 是已知的，且答案为 \(f_{2^n-1}\) 在 \(f\) 后半段中，于是就在 \(O(n^2 2^n)\) 完成了此题。

顺带一提，子集卷积求逆直接分别对 \(G_0\) , \(G_1\) FMT，对位做多项式除法，最后 IFMT 即可。由于项数很少只有 \(n\)，暴力除即可。

Code:

#include<bits/stdc++.h>
#define R register int
#define ll long long
#define I inline
using namespace std;
const int N=20,M=1<<N,P=1e9+7;
int n,l,c[N][N],g[M],pc[M];//pc即popcount，二进制中1的个数
I void pls(int &a,int b){a+=b;if(a>=P)a-=P;}
I void mns(int &a,int b){a-=b;if(a<0)a+=P;}
struct poly
{
	int a[N];
	void operator/=(const poly &p)
	{
		for(R i=1;i<=n;i++)
		{
			for(R j=0;j<i;j++)
				a[i]=(a[i]-1ll*a[j]*p.a[i-j])%P;
			if(a[i]<0)a[i]+=P;
		}
	}
	void operator+=(const poly &p){for(R i=0;i<=n;i++)pls(a[i],p.a[i]);}
	void operator-=(const poly &p){for(R i=0;i<=n;i++)mns(a[i],p.a[i]);}
};
poly g0[M>>1],g1[M>>1];
void FMT(poly *s)
{
	for(R i=1;i<l;i<<=1)
		for(R j=0;j<l;j+=i<<1)
			for(R k=j;k<j+i;k++)
				s[k+i]+=s[k];
}
void IFMT(poly *s)
{
	for(R i=1;i<l;i<<=1)
		for(R j=0;j<l;j+=i<<1)
			for(R k=j;k<j+i;k++)
				s[k+i]-=s[k];
}
int main()
{
	int rn,rl;
	scanf("%d",&rn);n=rn-1;//g, f拆成两半的长度处理
	l=1<<n;rl=1<<rn;
	for(R i=0;i<rn;i++)
		for(R j=0;j<rn;j++)
			scanf("%d",c[i]+j);
	g[0]=1;
	for(R i=1,j,k;i<rl;i++)
	{
		for(j=0;!(i>>j&1);j++);
		g[i]=g[i&i-1];
		for(k=j+1;i>>k;k++)
			if(i>>k&1)g[i]=g[i]*(c[j][k]+1ll)%P;
	}
	pc[0]=-1;//这里图方便把初值设为了-1，但其实后面会变成0
	for(R i=0;i<l;i++)
	{
		pc[i]=pc[i&i-1]+1;
		g0[i].a[pc[i]]=g[i];
		g1[i].a[pc[i]]=g[i|l];
		//注意这里占位多项式的处理，由于把g拆成了两段，后半段popcount应为原来的-1
	}
	FMT(g0);FMT(g1);
	for(R i=0;i<l;i++)g1[i]/=g0[i];
	IFMT(g1);
	printf("%d",g1[l-1].a[n]);
	return 0;
}