【JZOJ3303】城市规划

这篇博客讨论了如何解决一个城市规划问题,即如何在限制条件下确定无向连通图的不同构建方案数量。作者通过分析模数提出了分治与快速傅里叶变换(NTT)的方法,并详细解释了CDQ分治策略如何应用于求解这个问题,最终给出O(nlog^2n)的时间复杂度解决方案。

description

刚刚解决完电力网络的问题, 阿狸又被领导的任务给难住了.

刚才说过, 阿狸的国家有n 个城市, 现在国家需要在某些城市对之间建立一些贸易路线, 使得整个国家的任意两个城市都直接或间接的连通.

为了省钱, 每两个城市之间最多只能有一条直接的贸易路径. 对于两个建立路线的方案, 如果存在一个城市对, 在两个方案中是否建立路线不一样, 那么这两个方案就是不同的, 否则就是相同的. 现在你需要求出一共有多少不同的方案.

好了, 这就是困扰阿狸的问题. 换句话说, 你需要求出n 个点的简单(无重边无自环)无向连通图数目.

由于这个数字可能非常大, 你只需要输出方案数mod 1004535809(479 * 2 ^21 + 1)即可.


analysis

  • 看到模数就应该知道是分治 + N T T +NTT +NTT

  • f [ i ] 表 示 f[i]表示 f[i] i i i的答案,那么肯定是全部的方案数减去不合法的方案数

  • g [ i ] = 2 ( 2 i ) g[i]=2^{\left( \begin{matrix} 2\\ i\\ \end{matrix} \right)} g[i]=2(2i),也就是一共有 ( 2 i ) \left( \begin{matrix} 2\\ i\\ \end{matrix} \right) (2i)条边,就是 2 ( 2 i ) 2^{\left( \begin{matrix} 2\\ i\\ \end{matrix} \right)} 2(2i)种总方案数

  • 推推可以知道 f [ n ] = g [ n ] − ∑ i = 1 n − 1 ( i − 1 n − 1 ) g [ i ] ∗ f [ n − i ] f[n]=g[n]-\sum_{i=1}^{n-1}\left( \begin{matrix} i-1\\ n-1\\ \end{matrix} \right)g[i]*f[n-i] f[n]=g[n]i=1n1(i1n1)g[i]f[ni],拆组合数之后就是

  • f [ n ] = g [ n ] − ∑ i = 1 n − 1 ( n − 1 ) ! ( i − 1 ) ! ( n − i ) ! ∗ g [ i ] ∗ f [ n − i ] f[n]=g[n]-\sum_{i=1}^{n-1}{{(n-1)!}\over{(i-1)!(n-i)!}}*g[i]*f[n-i] f[n]=g[n]i=1n1(i1)!(ni)!(n1)!g[i]f[ni],移一下项就是

  • f [ n ] = g [ n ] − ( n − 1 ) ! ∑ i = 1 n − 1 g [ i ] ( i − 1 ) ! ∗ f [ n − i ] ( n − i ) ! f[n]=g[n]-(n-1)!\sum_{i=1}^{n-1}{{g[i]}\over{(i-1)!}}*{{f[n-i]}\over{(n-i)!}} f[n]=g[n](n1)!i=1n1(i1)!g[i](ni)!f[ni]

  • 后面那个就是卷积了,好像可以直接套 N T T NTT NTT,然而正着推是 O ( n 2 log ⁡ 2 n ) O(n^2\log_2n) O(n2log2n),我也不会多项式求逆

  • 于是考虑 C D Q CDQ CDQ分治,对于 [ l , r ] [l,r] [l,r]区间,单独求出所有 [ l , m i d ] [l,mid] [l,mid] [ m i d + 1 , r ] [mid+1,r] [mid+1,r]的贡献

  • 具体就是,把已经知道具体值的 f [ l . . m i d ] f[l..mid] f[l..mid]拉出来,与 g [ 1.. r − l + 1 ] g[1..r-l+1] g[1..rl+1]乘起来

  • 然后两个乘起来,后者长度是前者的两倍,多出来的那一段即为分别对 [ m i d + 1 , r ] [mid+1,r] [mid+1,r]每一位的贡献

  • 如果新多项式为 h h h h h h的第 i i i项系数为 f , g f,g f,g下标之和等于 i i i的系数积之和

  • 即比如 f [ 3 ] = g [ 3 ] − ( g [ 1 ] 0 ! ∗ f [ 2 ] 2 ! + g [ 2 ] 1 ! ∗ f [ 1 ] 1 ! ) f[3]=g[3]-({{g[1]}\over{0!}}*{{f[2]}\over{2!}}+{{g[2]}\over{1!}}*{{f[1]}\over{1!}}) f[3]=g[3](0!g[1]2!f[2]+1!g[2]1!f[1]),由于 f [ 2 ] f[2] f[2]已从 f [ 1 ] f[1] f[1]推得,于是

  • 注意为了方便下面我都把 f , g f,g f,g和阶乘的逆元并进 f ‘ , g ‘ f`,g` f,g里面了,不要看漏

  • ( f ‘ [ 1 ] ,   f ‘ [ 2 ] ,     0 ,     0 g ‘ [ 1 ] , g ‘ [ 2 ] , g ‘ [ 3 ] , g ‘ [ 4 ] h [ 1 ] , h [ 2 ] , h [ 3 ] , h [ 4 ] ) \left( \begin{matrix} f`[1],\ f`[2],\ \ \ 0,\ \ \ 0\\ g`[1],g`[2],g`[3],g`[4]\\ h[1],h[2],h[3],h[4]\\ \end{matrix} \right) f[1], f[2],   0,   0g[1],g[2],g[3],g[4]h[1],h[2],h[3],h[4],那么 h [ 3 ] = f ‘ [ 1 ] ∗ g ‘ [ 2 ] + f ‘ [ 2 ] ∗ g ‘ [ 1 ] h[3]=f`[1]*g`[2]+f`[2]*g`[1] h[3]=f[1]g[2]+f[2]g[1],恰好等于 f [ 3 ] − g [ 3 ] f[3]-g[3] f[3]g[3],第三位的答案就求出了

  • h [ 4 ] = f ‘ [ 1 ] ∗ g ‘ [ 3 ] + f ‘ [ 2 ] ∗ g ‘ [ 2 ] + 0 ∗ g ‘ [ 3 ] h[4]=f`[1]*g`[3]+f`[2]*g`[2]+0*g`[3] h[4]=f[1]g[3]+f[2]g[2]+0g[3] [ 1 , 2 ] [1,2] [1,2]区间对第四位的贡献就为 h [ 4 ] h[4] h[4]

  • 分治到 [ 3 , 3 ] [3,3] [3,3]区间时,这种方法同理会给第四位 f ‘ [ 3 ] ∗ g ‘ [ 1 ] f`[3]*g`[1] f[3]g[1]贡献,于是第四位的答案也求出了

  • 这样分治下去,可以知道是对的

  • 当区间长度为 1 1 1时,当前点前面的贡献都已经算完, f f f值即为 g [ l ] − ( l − 1 ) ! ∗ 贡 献 g[l]-(l-1)!*贡献 g[l](l1)!

  • 分治共 log ⁡ 2 n \log_2n log2n层,每层搞 N T T NTT NTT,复杂度即 O ( n log ⁡ 2 2 n ) O(n\log_2^2n) O(nlog22n)

  • 具体很难理解,可以多加思考,或者参考一下这里其实是我晚上做梦的时候想懂的

  • 不过其实多项式这块我还是超蒟,还是要多做题……


code

#include<stdio.h>
#include<string.h>
#include<algorithm>
#include<math.h>
#define g 3
#define mod 1004535809
#define MAXN 130000
#define MAXLEN MAXN*4
#define ll long long
#define fo(i,a,b) for (ll i=a;i<=b;++i)
#define fd(i,a,b) for (ll i=a;i>=b;--i)

using namespace std;

ll inv[MAXN+5],fac[MAXN+5],C[MAXLEN+5];
ll a[MAXLEN+5],b[MAXLEN+5],f[MAXLEN+5],rev[MAXLEN+5];
ll n,m,len,ans;

inline ll read()
{
	ll x=0,f=1;char ch=getchar();
	while (ch<'0' || '9'<ch){if (ch=='-')f=-1;ch=getchar();}
	while ('0'<=ch && ch<='9')x=x*10+ch-'0',ch=getchar();
	return x*f;
}
inline ll pow(ll x,ll y)
{
	ll z=1;
	while (y)
	{
		if (y&1)z=z*x%mod;
		y>>=1,x=x*x%mod;
	}
	return z;
}
inline ll work(ll x)
{
	ll y=1;
	while (y<x)y*=2;
	return y*2;
}
inline void init()
{
	fac[0]=fac[1]=1;
	fo(i,2,MAXN)fac[i]=fac[i-1]*i%mod;
	inv[MAXN]=pow(fac[MAXN],(ll)(mod-2))%mod;
	fd(i,MAXN,1)inv[i-1]=inv[i]*i%mod,C[i]=pow(2,i*(i-1)/2);
}
inline void ntt(ll a[],ll len,ll inv)
{
	ll bit=0;
	while ((1<<bit)<len)++bit;
	fo(i,0,len-1)
	{
		rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
		if (i<rev[i])swap(a[i],a[rev[i]]);
	}
	for (ll mid=1;mid<len;mid*=2)
	{
		ll tmp=pow(g,(mod-1)/(mid*2));
		if (inv==-1)tmp=pow(tmp,mod-2);
		for (ll i=0;i<len;i+=mid*2)
		{
			ll omega=1;
			for (ll j=0;j<mid;++j,omega=omega*tmp%mod)
			{
				ll x=a[i+j],y=omega*a[i+j+mid]%mod;
				a[i+j]=(x+y)%mod,a[i+j+mid]=(x-y+mod)%mod;
			}
		}
	}
}
inline void binary_search(ll l,ll r)
{
	if (l==r)
	{
		f[l]=(C[l]-f[l]*fac[l-1]%mod+mod)%mod;
		return;
	}
	ll mid=(l+r)>>1;
	binary_search(l,mid);
	ll len=work(r-l+1),invv=pow(len,mod-2);
	fo(i,0,len-1)a[i]=b[i]=0;
	fo(i,1,mid-l+1)a[i]=f[l+i-1]*inv[l+i-2]%mod;
	fo(i,1,r-l)b[i]=C[i]*inv[i]%mod;
	ntt(a,len,1),ntt(b,len,1);
	fo(i,0,len-1)a[i]=(a[i]*b[i])%mod;
	ntt(a,len,-1);
	fo(i,0,len-1)a[i]=(a[i]*invv)%mod;
	fo(i,mid+1,r)(f[i]+=a[i-l+1])%=mod;
	binary_search(mid+1,r);	
}
int main()
{
	n=read(),init(),m=work(n);
	binary_search(1,m/2);
	printf("%lld\n",f[n]);
	return 0;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值