我noip初赛最后半个小时就在数4个点的连通图有几个,结果还数错了…然后我就很膜拜这题,终于被我在jzoj找到了。
题意:n个不同的点的构成的连通图有多少种,模一个费马素数。n≤2e5
设f[i]为i个点的连通图方案数,思路为所有方案减去不合法的。
所有方案数为
2i∗(i−1)/2
。
对于不连通的,枚举1号点所在连通块大小为j。选这剩下j-1个点的方案数为
Cj−1i−1
。连通块内方案数为f[j],连通块外方案数为
2(i−j)∗(i−j−1)/2
。
整理一下,有
f[i]=2i∗(i−1)/2−∑j=1i−1f[j]∗Cj−1i−1∗2(i−j)∗(i−j−1)/2
两边同除(i-1)!,有
f[i](i−1)!=2i∗(i−1)/2(i−1)!−∑j=1i−1f[j]∗2(i−j)∗(i−j−1)/2(j−1)!∗(i−j)!
想到若1号点所在的连通块大小为i的含义,尝试将和式的j取i,发现它正好等于左边,有
∑j=1if[j]∗2(i−j)∗(i−j−1)/2(j−1)!∗(i−j)!=2i∗(i−1)/2(i−1)!
继续
∑j=1if[j](j−1)!∗2(i−j)∗(i−j−1)/2(i−j)!=2i∗(i−1)/2(i−1)!
左边已经是个卷积了。我们要求的f在左边,求个多项式逆元就可以了。
老是T的NTT
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define mmst(a, b) memset(a,b,sizeof(a))
#define mmcp(a, b) memcpy(a,b,sizeof(a))
using namespace std;
typedef long long LL;
const int N=400400;
const LL p=1004535809,g=3;
LL cheng(LL a,LL b)
{
LL res=1;
for(;b;b>>=1,a=a*a%p)
if(b&1)
res=res*a%p;
return res;
}
int n,rev[N];
void init(int lim)
{
n=1;
int k=-1;
while(n<lim)
n<<=1,k++;
for(int i=0;i<n;i++)
rev[i]=(rev[i>>1]>>1) | ((i&1)<<k);
}
void ntt(LL *a,bool ops)
{
for(int i=0;i<n;i++)
if(i<rev[i])
swap(a[i],a[rev[i]]);
for(int l=2;l<=n;l<<=1)
{
int m=l>>1;
LL wn;
if(ops)
wn=cheng(g,(p-1)/l);
else
wn=cheng(g,p-1-(p-1)/l);
for(int i=0;i<n;i+=l)
{
LL w=1;
for(int k=0;k<m;k++)
{
LL t=a[i+k+m]*w%p;
a[i+k+m]=(a[i+k]-t+p)%p;
a[i+k]=(a[i+k]+t)%p;
w=w*wn%p;
}
}
}
if(!ops)
{
LL Inv=cheng(n,p-2);
for(int i=0;i<n;i++)
a[i]=a[i]*Inv%p;
}
}
int nn;
LL jc[N],Ijc[N],B[N],C[N],D[N],E[N];
void ny(int x)
{
if(x==1)
return;
ny(x/2);
for(int i=0;i<x;i++)
E[i]=B[i];
init(x*2);
ntt(E,1);
ntt(D,1);
for(int i=0;i<n;i++)
D[i]=D[i]*(2-E[i]*D[i]%p+p)%p;
ntt(D,0);
for(int i=0;i<n;i++)
E[i]=0;
for(int i=x;i<n;i++)
D[i]=0;
}
int main()
{
cin>>nn;
jc[0]=1;
for(int i=1;i<=nn;i++)
jc[i]=jc[i-1]*i%p;
for(int i=0;i<=nn;i++)
Ijc[i]=cheng(jc[i],p-2);
for(LL i=1;i<=nn;i++)
C[i]=cheng(2,i*(i-1)/2)*Ijc[i-1]%p;
for(LL i=0;i<=nn;i++)
B[i]=cheng(2,i*(i-1)/2)*Ijc[i]%p;
B[0]=1;
D[0]=1;
init(nn+2);
ny(n);
init(nn*2+2);
ntt(D,1);
ntt(C,1);
for(int i=0;i<n;i++)
C[i]=C[i]*D[i]%p;
ntt(C,0);
cout<<C[nn]*jc[nn-1]%p<<endl;
return 0;
}