指数型生成函数&[bzoj3456]城市规划

11 篇文章 0 订阅
1 篇文章 0 订阅

前言

打完多项式板子后的第一题+清真的题意

题目相关

题目链接

题目大意

n n n个点的简单(无重边无自环)无向连通图数目
输出模 1004535809 ( 479 ∗ 2 21 + 1 ) 1004535809(479*2^{21}+1) 1004535809(479221+1)

数据范围

n ≤ 130000 n\le 130000 n130000

题解

暴力

我只推出了个暴力,感觉和多项式没啥关系
如果没有任何限制,那么答案为 2 C n 2 2^{C_n^2} 2Cn2
F n F_n Fn为答案,写出递推式(用所有的方案数减去不连通的方案数,枚举最后一个点所处的联通块大小即可)
F i = 2 C i 2 − ∑ j = 1 i − 1 F j ∗ C i − 1 j − 1 ∗ 2 C i − j 2 F_i=2^{C_i^2}-\sum_{j=1}^{i-1}F_j*C_{i-1}^{j-1}*2^{C_{i-j}^2} Fi=2Ci2j=1i1FjCi1j12Cij2
递推+快速幂
复杂度 O ( n 2 l o g p ) \mathcal O(n^2logp) O(n2logp)
要是没有代码长度限制就可以打表过了(大雾)

正解

和暴力几乎没啥关系
考虑指数型生成函数
我们来观察指数型生成函数的性质,列出两个指数型生成函数
A ( x ) = ∑ i = 0 ∞ a i x i i ! A(x)=\sum_{i=0}^{\infty}a_i\frac{x^i}{i!} A(x)=i=0aii!xi
B ( x ) = ∑ i = 0 ∞ b i x i i ! B(x)=\sum_{i=0}^{\infty}b_i\frac{x^i}{i!} B(x)=i=0bii!xi
我们观察一下 A ( x ) ⋅ B ( x ) A(x)·B(x) A(x)B(x)
A ( x ) ⋅ B ( X ) = ( ∑ i = 0 ∞ a i x i i ! ) ⋅ ( ∑ i = 0 ∞ b i x i i ! ) = ∑ i = 0 ∞ a i x i i ! ∑ j = 0 ∞ b j x j j ! = ∑ i = 0 ∞ ∑ j = 0 ∞ a i x i i ! b j x j j ! = ∑ n = 0 ∞ ( ∑ i = 0 ∞ ( n i ) a i b n − i ) x n n ! \begin{aligned} A(x)·B(X)&=(\sum_{i=0}^{\infty}a_i\frac{x^i}{i!})·(\sum_{i=0}^{\infty}b_i\frac{x^i}{i!})\\ &=\sum_{i=0}^{\infty}a_i\frac{x^i}{i!}\sum_{j=0}^{\infty}b_j\frac{x^j}{j!}\\ &=\sum_{i=0}^{\infty}\sum_{j=0}^{\infty}a_i\frac{x^i}{i!}b_j\frac{x^j}{j!}\\ &=\sum_{n=0}^{\infty}(\sum_{i=0}^{\infty}\binom{n}{i}a_ib_{n-i})\frac{x^n}{n!}\\ \end{aligned} A(x)B(X)=(i=0aii!xi)(i=0bii!xi)=i=0aii!xij=0bjj!xj=i=0j=0aii!xibjj!xj=n=0(i=0(in)aibni)n!xn
发现这个形式非常的奇妙
如果 a i a_i ai代表在第一个集合中选 i i i个的方案数, b i b_i bi代表在第二个集合中选 i i i个的方案数
那么相乘后我们就可以知道总共选 i i i个的方案数

考虑划分
我们发现,如果我们按照划分出一个大小为 i i i的集合的方案数建一个指数型生成函数 F ( x ) F(x) F(x)
那么我们就能计算出对于 n n n个元素划分成 m m m个集合的方案数
对于 m = 0 m=0 m=0,答案的多项式为 1 1 1(即只有 n = 0 n=0 n=0时方案数为 1 1 1
对于 m = 1 m=1 m=1,答案多项式为 F ( x ) 1 ! \frac{F(x)}{1!} 1!F(x)
对于 m = 2 m=2 m=2,答案多项式为 F 2 ( x ) 2 ! \frac{F^2(x)}{2!} 2!F2(x)
对于 m = 3 m=3 m=3,答案多项式为 F 3 ( x ) 3 ! \frac{F^3(x)}{3!} 3!F3(x)
···
由于划分先后的问题,所以要除以阶乘
如果我们要求每个 n n n的划分总方案数
F 0 ( x ) 0 ! + F 1 ( x ) 1 ! + F ( x ) 2 ! ⋅ ⋅ ⋅ = ∑ i = 0 ∞ F i ( x ) i ! \frac{F^0(x)}{0!}+\frac{F^1(x)}{1!}+\frac{F^(x)}{2!}···=\sum_{i=0}^\infty \frac{F^i(x)}{i!} 0!F0(x)+1!F1(x)+2!F(x)=i=0i!Fi(x)
根据底数为 e e e的泰勒展开式 e x = ∑ i = 0 ∞ x i i ! e^x=\sum_{i=0}^\infty\frac {x^i}{i!} ex=i=0i!xi
我们得到 ∑ i = 0 ∞ F i ( x ) i ! = e F ( x ) \sum_{i=0}^\infty \frac{F^i(x)}{i!}=e^{F(x)} i=0i!Fi(x)=eF(x)


回到本题: 我们思考,在这一题中,所有联通图都是无限制图的一个划分
F ( x ) F(x) F(x)为联通图方案数的多项式, G ( x ) G(x) G(x)为无限制图方案数的多项式
容易发现 e F ( x ) = G ( x ) e^{F(x)}=G(x) eF(x)=G(x)
所以 F ( x ) = l n ( G ( x ) ) F(x)=ln(G(x)) F(x)=ln(G(x))
我们发现 G ( x ) G(x) G(x)非常好求(暴力的那部分已经说过了)
然后直接多项式ln即可

代码

拉了板子所以代码很长

#include<cstdio>
#include<cctype>
#include<cstring>
#include<cstdlib>
#include<vector>
namespace fast_IO
{
    const int IN_LEN=10000000,OUT_LEN=10000000;
    char ibuf[IN_LEN],obuf[OUT_LEN],*ih=ibuf+IN_LEN,*oh=obuf,*lastin=ibuf+IN_LEN,*lastout=obuf+OUT_LEN-1;
    inline char getchar_(){return (ih==lastin)&&(lastin=(ih=ibuf)+fread(ibuf,1,IN_LEN,stdin),ih==lastin)?EOF:*ih++;}
    inline void putchar_(const char x){if(oh==lastout)fwrite(obuf,1,oh-obuf,stdout),oh=obuf;*oh++=x;}
    inline void flush(){fwrite(obuf,1,oh-obuf,stdout);}
}
using namespace fast_IO;
#define getchar() getchar_()
#define putchar(x) putchar_((x))
typedef long long ll;
#define rg register
template <typename T> inline T max(const T a,const T b){return a>b?a:b;}
template <typename T> inline T min(const T a,const T b){return a<b?a:b;}
template <typename T> inline T mind(T&a,const T b){a=a<b?a:b;}
template <typename T> inline T maxd(T&a,const T b){a=a>b?a:b;}
template <typename T> inline T abs(const T a){return a>0?a:-a;}
template <typename T> inline void swap(T&a,T&b){T c=a;a=b;b=c;}
template <typename T> inline void swap(T*a,T*b){T c=a;a=b;b=c;}
template <typename T> inline T gcd(const T a,const T b){if(!b)return a;return gcd(b,a%b);}
template <typename T> inline T square(const T x){return x*x;};
template <typename T> inline void read(T&x)
{
    char cu=getchar();x=0;bool fla=0;
    while(!isdigit(cu)){if(cu=='-')fla=1;cu=getchar();}
    while(isdigit(cu))x=x*10+cu-'0',cu=getchar();
    if(fla)x=-x;  
}
template <typename T> void printe(const T x)
{
    if(x>=10)printe(x/10);
    putchar(x%10+'0');
}
template <typename T> inline void print(const T x)
{
    if(x<0)putchar('-'),printe(-x);
    else printe(x);
}
const int maxn=2097152,mod=1004535809;
inline int Md(const int x){return x>=mod?x-mod:x;}
template<typename T>
inline int pow(int x,T y)
{
    rg int res=1;x%=mod;
    for(;y;y>>=1,x=(ll)x*x%mod)if(y&1)res=(ll)res*x%mod;
    return res;
}
namespace Poly///namespace of Poly
{
int W_[maxn],FW_[maxn],ha[maxn],hb[maxn],Inv[maxn];
inline void init(const int x)
{
    rg int tim=0,lenth=1;
    while(lenth<x)lenth<<=1,tim++;
    for(rg int i=1;i<lenth;i++)
    {
        W_[i]=pow(3,(mod-1)/i/2);
        FW_[i]=pow(W_[i],mod-2);
    }
    Inv[0]=Inv[1]=1;
    for(rg int i=2;i<x;i++)Inv[i]=(ll)(mod-mod/i)*Inv[mod%i]%mod;
}
int L;
inline void NTT(int*A,const int fla)//prepare:init L 
{
    for(rg int i=0,j=0;i<L;i++)
    {
        if(i>j)swap(A[i],A[j]);
        for(rg int k=L>>1;(j^=k)<k;k>>=1);
    }
    for(rg int i=1;i<L;i<<=1)
    {
        const int w=fla==-1?FW_[i]:W_[i];
        for(rg int j=0,J=i<<1;j<L;j+=J)
        {
            int K=1;
            for(rg int k=0;k<i;k++,K=(ll)K*w%mod)
            {
                const int x=A[j+k],y=(ll)A[j+k+i]*K%mod;
                A[j+k]=Md(x+y),A[j+k+i]=Md(mod+x-y);
            }
        }
    }
}
inline int Quadratic_residue(const int a)
{
	if(a==0)return 0;
	int b=(rand()<<14^rand())%mod;
	while(pow(b,(mod-1)>>1)!=mod-1)b=(rand()<<14^rand())%mod;
	int s=mod-1,t=0,x,inv=pow(a,mod-2),f=1;
	while(!(s&1))s>>=1,t++,f<<=1;
	t--,x=pow(a,(s+1)>>1),f>>=1;
	while(t)
	{
		f>>=1;
		if(pow((int)((ll)inv*x%mod*x%mod),f)!=1)x=(ll)x*pow(b,s)%mod;
		t--,s<<=1;
	}
	return min(x,mod-x);
}
struct poly
{
    std::vector<int>A;
    poly(){A.resize(0);}
    poly(const int x){A.resize(1),A[0]=x;}
    inline int&operator[](const int x){return A[x];}
    inline int operator[](const int x)const{return A[x];}
    inline void clear(){A.clear();}
    inline unsigned int size()const{return A.size();}
    inline void resize(const unsigned int x){A.resize(x);}
    void RE(const int x)
    {
        A.resize(x);
        for(rg int i=0;i<x;i++)A[i]=0; 
    }
    void readin(const int MAX)
    {
        A.resize(MAX);
        for(rg int i=0;i<MAX;i++)read(A[i]);
    }
    void putout()const
    {
        for(rg unsigned int i=0;i<A.size();i++)print(A[i]),putchar(' ');
    }
    inline poly operator +(const poly b)const
    {
        poly RES;
        RES.resize(max(size(),b.size()));
        for(rg unsigned int i=0;i<RES.size();i++)RES[i]=Md((i<size()?A[i]:0)+(i<b.size()?b[i]:0));
        return RES;
    }
    inline poly operator -(const poly b)const
    {
        poly RES;
        RES.resize(max(size(),b.size()));
        for(rg unsigned int i=0;i<RES.size();i++)RES[i]=Md((i<size()?A[i]:0)+mod-(i<b.size()?b[i]:0));
        return RES;
    }
    inline poly operator *(const int b)const
	{
		poly RES=*this;
        for(rg unsigned int i=0;i<RES.size();i++)RES[i]=(ll)RES[i]*b%mod;
        return RES;
	}
	inline poly operator /(const int b)const
	{
		poly RES=(*this)*pow(b,mod-2);
    	return RES;
	}
    inline poly operator *(const poly b)const
    {
        const int RES=A.size()+b.size()+1;
        L=1;while(L<RES)L<<=1;
        poly c;c.A.resize(RES);
        memset(ha,0,sizeof(int)*L);
        memset(hb,0,sizeof(int)*L);
        for(rg unsigned int i=0;i<A.size();i++)ha[i]=A[i];
        for(rg unsigned int i=0;i<b.A.size();i++)hb[i]=b.A[i];
        NTT(ha,1),NTT(hb,1);
        for(rg int i=0;i<L;i++)ha[i]=(ll)ha[i]*hb[i]%mod;
        NTT(ha,-1);
        const int inv=pow(L,mod-2);
        for(rg int i=0;i<RES;i++)c.A[i]=(ll)ha[i]*inv%mod;
        return c;
    }
    inline poly inv()const
    {
        poly C;
        if(A.size()==1){C=*this;C[0]=pow(C[0],mod-2);return C;}
        C.resize((A.size()+1)>>1);
        for(rg unsigned int i=0;i<C.size();i++)C[i]=A[i];
        C=C.inv();
        L=1;while(L<(int)size()*2-1)L<<=1;
        for(rg unsigned int i=0;i<A.size();i++)ha[i]=A[i];
        for(rg unsigned int i=0;i<C.size();i++)hb[i]=C[i];
        memset(ha+A.size(),0,sizeof(int)*(L-A.size()));
        memset(hb+C.size(),0,sizeof(int)*(L-C.size()));
        NTT(ha,1),NTT(hb,1);
        for(rg int i=0;i<L;i++)ha[i]=(ll)(2-(ll)hb[i]*ha[i]%mod+mod)*hb[i]%mod;
        NTT(ha,-1);
        const int inv=pow(L,mod-2);
        C.resize(size());
        for(rg unsigned int i=0;i<size();i++)C[i]=(ll)ha[i]*inv%mod;
        return C;
    }
/*    inline poly inv()const
    {
        poly C;
        if(A.size()==1){C=*this;C[0]=pow(C[0],mod-2);return C;}
        C.resize((A.size()+1)>>1);
        for(rg unsigned int i=0;i<C.size();i++)C[i]=A[i];
        C=C.inv();
        poly D=(poly)2-*this*C;
        D.resize(size());
        D=D*C;
        D.resize(size());
        return D;
    }*///大常数版本 
    inline void Reverse(const int n)
    {
    	A.resize(n);
    	for(rg int i=0,j=n-1;i<j;i++,j--)swap(A[i],A[j]);
    }
    inline poly operator /(const poly B)const
	{
		poly a=*this,b=B;a.Reverse(size()),b.Reverse(B.size());
		b.resize(size()-B.size()+1);
		b=b.inv();
		b=b*a;
		b.Reverse(size()-B.size()+1);
		return b;
	}
	inline poly operator %(const poly B)const
	{
		poly C=(*this)-(*this)/B*B;C.resize(B.size()-1);
		return C;
	}
	inline poly diff()const
	{
		poly C;C.resize(size()-1);
		for(rg unsigned int i=1;i<size();i++)C[i-1]=(ll)A[i]*i%mod;
		return C;
	}
	inline poly inte()const
	{
		poly C;C.resize(size()+1);
		for(rg unsigned int i=0;i<size();i++)C[i+1]=(ll)A[i]*Inv[i+1]%mod;
		C[0]=0;
		return C;
	}
	inline poly ln()const
	{
		poly C=(diff()*inv()).inte();C.resize(size());
		return C;
	}
	inline poly exp()const
	{
		poly C;
		if(size()==1){C=*this;C[0]=1;return C;}
		C.resize((size()+1)>>1);
		for(rg unsigned int i=0;i<C.size();i++)C[i]=A[i];
		C=C.exp();C.resize(size());
		poly D=(poly)1-C.ln()+*this;
		D=D*C;
		D.resize(size());
		return D;
	}
	inline poly sqrt()const
	{
		poly C;
		if(size()==1)
		{
			C=*this;C[0]=Quadratic_residue(C[0]);
			return C;
		}
		C.resize((size()+1)>>1);
		for(rg unsigned int i=0;i<C.size();i++)C[i]=A[i];
		C=C.sqrt();C.resize(size());
		C=(C+*this*C.inv())*(int)499122177;
		C.resize(size());
		return C;
	}
	inline poly operator >>(const unsigned int x)const
    {
    	poly C;if(size()<x){C.resize(0);return C;}
		C.resize(size()-x);
    	for(rg unsigned int i=0;i<C.size();i++)C[i]=A[i+x];
    	return C;
    }
    inline poly operator <<(const unsigned int x)const
    {
    	poly C;C.RE(size()+x);
    	for(rg unsigned int i=0;i<size();i++)C[i+x]=A[i];
    	return C;
    }
    inline poly Pow(const unsigned int x)const
    {
    	for(rg unsigned int i=0;i<size();i++)
			if(A[i])
			{
				poly C=((((*this/A[i])>>i).ln()*x).exp()*pow(A[i],x))<<(min(i*x,size()));
				C.resize(size());
				return C;
			}
    	return *this;
    }
};
}///namespace of Poly
Poly::poly a;
int n;
int fac[maxn+1],inv[maxn+1]; 
int main()
{
	fac[0]=1;for(rg int i=1;i<=maxn;i++)fac[i]=(ll)fac[i-1]*i%mod;
	inv[maxn]=pow(fac[maxn],mod-2);for(rg int i=maxn;i>=1;i--)inv[i-1]=(ll)inv[i]*i%mod;
    Poly::init(maxn);///namespace of Poly
    read(n);
    a.RE(n+1);
	for(rg int i=0;i<=n;i++)a[i]=(ll)pow(2,((ll)i*(i-1))>>1)*inv[i]%mod;
	a=a.ln();
	print((ll)a[n]*fac[n]%mod);
	return flush(),0;
}

总结

好推+好用的指数型生成函数!

  • 1
    点赞
  • 2
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值