zoukankan      html  css  js  c++  java
  • 【BZOJ3456】轩辕朗的城市规划 无向连通图计数 CDQ分治 FFT 多项式求逆 多项式ln

    题解

    分治FFT

      设(f_i)(i)个点组成的无向图个数,(g_i)(i)个点组成的无向连通图个数

      经过简单的推导(枚举(1)所在的连通块大小),有:

    [f_i=2^{frac{i(i-1)}{2}} ]

    [egin{align} g_i&=f_i-sum_{j=1}^{i-1}inom{n-1}{j-1}g_jf_{i-j}\ &=f_i-(i-1)!sum_{j=1}^{i-1}frac{g_j}{(j-1)!}frac{f_{i-j}}{(i-j)!} end{align} ]

      用CDQ分治+FFT优化。

      就是每次先求出(g_lcdots g_{mid+1}),然后卷上(f_1cdots f_{len}),加到(g_{mid+1}cdots g_r)上面去。

      时间复杂度:(O(nlog^2 n))

    多项式求逆

    [egin{align} f_i&=sum_{j=1}^ninom{i-1}{j-1}g_jf_{i-j}\ frac{f_i}{(i-1)!}&=sum_{j=1}^nfrac{g_j}{(j-1)!}frac{f_{i-j}}{(i-j)!} end{align} ]

      设

    [egin{align} A&=sum_{igeq 1}frac{f_i}{(i-1)!}x^i\ B&=sum_{igeq 1}frac{g_i}{(i-1)!}x^i\ C&=sum_{igeq 0}frac{f_i}{i!}x^i end{align} ]

      所以

    [egin{align} A&=B imes C~~~~~~(mod~x^{n+1})\ B&=A imes C^{-1}~~(mod~x^{n+1}) end{align} ]

      时间复杂度:(O(nlog n))

    多项式求ln

      设

    [egin{align} G(x)&=sum_{igeq 0}frac{2^{inom{i}{2}}}{i!}x^i\ F(x)&=sum_{igeq 0}frac{f_i}{i!}x^i end{align} ]

      根据指数生成函数和有标号计数的那套理论,
      由(0)个连通块组成的连通图的个数为(frac{{F(x)}^0}{0!})
      由(1)个连通块组成的连通图的个数为(frac{{F(x)}^1}{1!})
      由(2)个连通块组成的连通图的个数为(frac{{F(x)}^2}{2!})
      (vdots)

      这些加起来就是无向图的个数(G(x))

      所以

    [egin{align} G(x)&=frac{{F(x)}^0}{0!}+frac{{F(x)}^1}{1!}+frac{{F(x)}^2}{2!}+frac{{F(x)}^3}{3!}+cdots\ &=sum_{igeq 0}frac{{F(x)}^i}{i!}\ &=e^{F(x)}\ F(x)&=ln(G(x)) end{align} ]

      直接求ln即可。

      时间复杂度:(O(nlog n))

    代码

    分治FFT

    #include<cstdio>
    #include<cstring>
    #include<algorithm>
    #include<cstdlib>
    #include<ctime>
    #include<utility>
    using namespace std;
    typedef long long ll;
    typedef pair<int,int> pii;
    ll p=1004535809;
    ll fp(ll a,ll b)
    {
    	ll s=1;
    	while(b)
    	{
    		if(b&1)
    			s=s*a%p;
    		a=a*a%p;
    		b>>=1;
    	}
    	return s;
    }
    ll inv(ll x)
    {
    	return fp(x,p-2);
    }
    namespace ntt
    {
    	ll w1[300010];
    	ll w2[300010];
    	int rev[300010];
    	int n;
    	void init(int m)
    	{
    		n=1;
    		while(n<=m)
    			n<<=1;
    		int i;
    		for(i=2;i<=n;i<<=1)
    		{
    			w1[i]=fp(3,(p-1)/i);
    			w2[i]=inv(w1[i]);
    		}
    		rev[0]=0;
    		for(i=1;i<n;i++)
    			rev[i]=(rev[i>>1]>>1)|(i&1?n>>1:0);
    	}
    	void ntt(ll *a,int t)
    	{
    		int i,j,k;
    		ll u,v,w,wn;
    		for(i=0;i<n;i++)
    			if(rev[i]<i)
    				swap(a[i],a[rev[i]]);
    		for(i=2;i<=n;i<<=1)
    		{
    			wn=(t==1?w1[i]:w2[i]);
    			for(j=0;j<n;j+=i)
    			{
    				w=1;
    				for(k=j;k<j+i/2;k++)
    				{
    					u=a[k];
    					v=a[k+i/2]*w%p;
    					a[k]=(u+v)%p;
    					a[k+i/2]=(u-v+p)%p;
    					w=w*wn%p;
    				}
    			}
    		}
    		if(t==-1)
    		{
    			u=inv(n);
    			for(i=0;i<n;i++)
    				a[i]=a[i]*u%p;
    		}
    	}
    };
    ll a[300010];
    ll b[300010];
    ll fac[300010];
    ll invfac[300010];
    ll in[300010];
    ll f[300010];
    ll g[300010];
    void solve(int l,int r)
    {
    	if(l==r)
    	{
    		g[l]=(f[l]-fac[l-1]*g[l]%p+p)%p;
    		return;
    	}
    	int mid=(l+r)>>1;
    	solve(l,mid);
    	int len=r-l+1;
    	ntt::init(len);
    	int i;
    	a[0]=b[0]=0;
    	if(l==3&&r==4)
    		int x=1;
    	if(l==1&&r==4)
    		int x=1;
    	for(i=1;i<=mid-l+1;i++)
    		a[i]=g[i+l-1]*invfac[i+l-1-1]%p;
    	for(i=mid-l+1+1;i<ntt::n;i++)
    		a[i]=0;
    	for(i=1;i<=len;i++)
    		b[i]=f[i]*invfac[i]%p;
    	for(i=len+1;i<ntt::n;i++)
    		b[i]=0;
    	ntt::ntt(a,1);
    	ntt::ntt(b,1);
    	for(i=0;i<ntt::n;i++)
    		a[i]=a[i]*b[i]%p;
    	ntt::ntt(a,-1);
    	for(i=mid+1;i<=r;i++)
    		g[i]=(g[i]+a[i-l+1])%p;
    	solve(mid+1,r);
    //	fprintf(stderr,"%d %d
    ",l,r);
    }
    int main()
    {
    //	freopen("bzoj3456.in","r",stdin);
    //	freopen("bzoj3456.out","w",stdout);
    	int n;
    	scanf("%d",&n);
    	invfac[0]=fac[0]=1;
    	int i;
    	in[0]=in[1]=1;
    	for(i=2;i<=n;i++)
    		in[i]=(-(p/i)*in[p%i]%p+p)%p;
    	for(i=1;i<=n;i++)
    	{
    		fac[i]=fac[i-1]*i%p;
    		invfac[i]=invfac[i-1]*in[i]%p;
    	}
    	for(i=1;i<=n;i++)
    		f[i]=fp(2,(ll(i)*(i-1)/2)%(p-1));
    	memset(g,0,sizeof g);
    	solve(1,n);
    	printf("%lld
    ",g[n]);
    	return 0;
    }
    

    多项式求逆

    #include<cstdio>
    #include<cstring>
    #include<algorithm>
    #include<cstdlib>
    #include<ctime>
    #include<utility>
    using namespace std;
    typedef long long ll;
    typedef pair<int,int> pii;
    ll p=1004535809;
    ll g=3;
    ll fp(ll a,ll b)
    {
    	ll s=1;
    	while(b)
    	{
    		if(b&1)
    			s=s*a%p;
    		a=a*a%p;
    		b>>=1;
    	}
    	return s;
    }
    namespace ntt
    {
    	ll w1[300010];
    	ll w2[300010];
    	int rev[300010];
    	int n;
    	void init(int m)
    	{
    		n=1;
    		while(n<m)
    			n<<=1;
    		int i;
    		for(i=2;i<=n;i++)
    		{
    			w1[i]=fp(g,(p-1)/i);
    			w2[i]=fp(w1[i],p-2);
    		}
    		rev[0]=0;
    		for(i=1;i<n;i++)
    			rev[i]=(rev[i>>1]>>1)|(i&1?n>>1:0);
    	}
    	void ntt(ll *a,int t)
    	{
    		int i,j,k;
    		ll u,v,w,wn;
    		for(i=0;i<n;i++)
    			if(rev[i]<i)
    				swap(a[i],a[rev[i]]);
    		for(i=2;i<=n;i<<=1)
    		{
    			wn=(t==1?w1[i]:w2[i]);
    			for(j=0;j<n;j+=i)
    			{
    				w=1;
    				for(k=j;k<j+i/2;k++)
    				{
    					u=a[k];
    					v=a[k+i/2]*w%p;
    					a[k]=(u+v)%p;
    					a[k+i/2]=(u-v+p)%p;
    					w=w*wn%p;
    				}
    			}
    		}
    		if(t==-1)
    		{
    			u=fp(n,p-2);	
    			for(i=0;i<n;i++)
    				a[i]=a[i]*u%p;
    		}
    	}
    	ll x[300010];
    	ll y[300010];
    	void copy_clear(ll *a,ll *b,int m)
    	{
    		int i;
    		for(i=0;i<m;i++)
    			a[i]=b[i];
    		for(i=m;i<n;i++)
    			a[i]=0;
    	}
    	void copy(ll *a,ll *b,int m)
    	{
    		int i;
    		for(i=0;i<m;i++)
    			a[i]=b[i];
    	}
    	void inverse(ll *a,ll *b,int m)
    	{
    		if(m==1)
    		{
    			b[0]=fp(a[0],p-2);
    			return;
    		}
    		inverse(a,b,m>>1);
    		init(2*m);
    		copy_clear(x,a,m);
    		copy_clear(y,b,m>>1);
    		ntt(x,1);
    		ntt(y,1);
    		int i;
    		for(i=0;i<n;i++)
    			x[i]=(2*y[i]%p-x[i]*y[i]%p*y[i]%p+p)%p;
    		ntt(x,-1);
    		copy(b,x,m);
    	}
    };
    ll a[300010];
    ll b[300010];
    ll c[300010];
    ll fac[300010];
    int main()
    {
    //	freopen("bzoj3456.in","r",stdin);
    	int n;
    	scanf("%d",&n);
    	int i;
    	fac[0]=1;
    	int m=1;
    	while(m<=2*n)
    		m<<=1;
    	for(i=1;i<=n;i++)
    		fac[i]=fac[i-1]*i%p;
    	a[0]=0;
    	for(i=1;i<=n;i++)
    		a[i]=fp(2,(ll(i-1)*i/2)%(p-1))*fp(fac[i-1],p-2)%p;
    	b[0]=1;
    	for(i=1;i<=n;i++)
    		b[i]=fp(2,(ll(i-1)*i/2)%(p-1))*fp(fac[i],p-2)%p;
    	ntt::inverse(b,c,m>>1);
    	for(i=n+1;i<m;i++)
    		c[i]=0;
    	ntt::init(m);
    	ntt::ntt(a,1);
    	ntt::ntt(c,1);
    	for(i=0;i<m;i++)
    		a[i]=a[i]*c[i]%p;
    	ntt::ntt(a,-1);
    	printf("%d
    ",a[n]*fac[n-1]%p);
    	return 0;
    }
    

    多项式ln

    #include<cstdio>
    #include<cstring>
    #include<algorithm>
    #include<cstdlib>
    #include<ctime>
    #include<utility>
    #include<cmath>
    #include<functional>
    using namespace std;
    typedef long long ll;
    typedef unsigned long long ull;
    typedef pair<int,int> pii;
    typedef pair<ll,ll> pll;
    void sort(int &a,int &b)
    {
    	if(a>b)
    		swap(a,b);
    }
    void open(const char *s)
    {
    #ifndef ONLINE_JUDGE
    	char str[100];
    	sprintf(str,"%s.in",s);
    	freopen(str,"r",stdin);
    	sprintf(str,"%s.out",s);
    	freopen(str,"w",stdout);
    #endif
    }
    int rd()
    {
    	int s=0,c;
    	while((c=getchar())<'0'||c>'9');
    	do
    	{
    		s=s*10+c-'0';
    	}
    	while((c=getchar())>='0'&&c<='9');
    	return s;
    }
    int upmin(int &a,int b)
    {
    	if(b<a)
    	{
    		a=b;
    		return 1;
    	}
    	return 0;
    }
    int upmax(int &a,int b)
    {
    	if(b>a)
    	{
    		a=b;
    		return 1;
    	}
    	return 0;
    }
    const ll p=1004535809;
    const ll g=3;
    ll fp(ll a,ll b)
    {
        ll s=1;
        while(b)
        {
            if(b&1)
                s=s*a%p;
            a=a*a%p;
            b>>=1;
        }
        return s;
    }
    const int maxn=600000;
    ll inv[maxn];
    namespace ntt
    {
        ll w1[maxn];
        ll w2[maxn];
        int rev[maxn];
        int n;
        void init(int m)
        {
            n=1;
            while(n<m)
                n<<=1;
            int i;
            for(i=2;i<=n;i++)
            {
                w1[i]=fp(g,(p-1)/i);
                w2[i]=fp(w1[i],p-2);
            }
            rev[0]=0;
            for(i=1;i<n;i++)
                rev[i]=(rev[i>>1]>>1)|(i&1?n>>1:0);
        }
        void ntt(ll *a,int t)
        {
            int i,j,k;
            ll u,v,w,wn;
            for(i=0;i<n;i++)
                if(rev[i]<i)
                    swap(a[i],a[rev[i]]);
            for(i=2;i<=n;i<<=1)
            {
                wn=(t==1?w1[i]:w2[i]);
                for(j=0;j<n;j+=i)
                {
                    w=1;
                    for(k=j;k<j+i/2;k++)
                    {
                        u=a[k];
                        v=a[k+i/2]*w%p;
    					a[k]=(u+v)%p;
    					a[k+i/2]=(u-v)%p;
                        w=w*wn%p;
                    }
                }
            }
            if(t==-1)
            {
                u=fp(n,p-2);    
                for(i=0;i<n;i++)
                    a[i]=a[i]*u%p;
            }
        }
        ll x[maxn];
        ll y[maxn];
        void copy_clear(ll *a,ll *b,int m)
        {
            int i;
            for(i=0;i<m;i++)
                a[i]=b[i];
            for(i=m;i<n;i++)
                a[i]=0;
        }
        void copy(ll *a,ll *b,int m)
        {
            int i;
            for(i=0;i<m;i++)
                a[i]=b[i];
        }
        void inverse(ll *a,ll *b,int m)
        {
            if(m==1)
            {
                b[0]=fp(a[0],p-2);
                return;
            }
            inverse(a,b,m>>1);
            init(m<<1);
            copy_clear(x,a,m);
            copy_clear(y,b,m>>1);
            ntt(x,1);
            ntt(y,1);
            int i;
            for(i=0;i<n;i++)
                x[i]=(2*y[i]-x[i]*y[i]%p*y[i])%p;
        	ntt(x,-1);
        	copy(b,x,m);
        }
        void integrate(ll *a,ll *b,int m)
    	{
    		int i;
    		for(i=0;i<m-1;i++)
    			b[i]=(i+1)*a[i+1]%p;
    		b[m-1]=0;
    	}
        void differential(ll *a,ll *b,int m)
        {
        	int i;
        	for(i=m-1;i>=1;i--)
        		b[i]=a[i-1]*inv[i]%p;
        	b[0]=0;
        }
        void ln(ll *a,ll *b,int m)
        {
        	static ll c[maxn],d[maxn],e[maxn];
        	integrate(a,c,m);
        	inverse(a,d,m);
        	init(m<<1);
        	ntt(c,1);
        	ntt(d,1);
        	int i;
        	for(i=0;i<n;i++)
        		e[i]=c[i]*d[i];
        	ntt(e,-1);
        	differential(e,b,m);
        }
        void exp(ll *a,ll *b,int m)
        {
        	if(m==1)
        	{
        		b[0]=1;
        		return;
        	}
        	exp(a,b,m>>1);
        	int i;
        	for(i=m>>1;i<m;i++)
        		b[i]=0;
        	ln(b,y,m);
        	init(m<<1);
        	copy_clear(x,a,m);
        	x[0]++;
        	for(i=0;i<m;i++)
        		x[i]=(x[i]-y[i])%p;
        	copy_clear(y,b,m);
        	ntt(x,1);
        	ntt(y,1);
        	for(i=0;i<n;i++)
        		x[i]=x[i]*y[i]%p;
        	ntt(x,-1);
        	for(i=0;i<m;i++)
        		b[i]=x[i];
        }
    };
    ll a[maxn];
    ll b[maxn];
    ll fac[maxn];
    ll invfac[maxn];
    int main()
    {
    	open("bzoj3456");
    	int i,n;
    	scanf("%d",&n);
    	int m=1;
    	while(m<=n+1)
    		m<<=1;
    	inv[0]=inv[1]=fac[0]=fac[1]=invfac[0]=invfac[1]=1;
    	for(i=2;i<m;i++)
    	{
    		inv[i]=-(p/i)*inv[p%i]%p;
    		fac[i]=fac[i-1]*i%p;
    		invfac[i]=invfac[i-1]*inv[i]%p;
    	}
    	for(i=0;i<m;i++)
    		a[i]=fp(2,(ll(i)*(i-1)/2)%(p-1))*invfac[i]%p;
    	ntt::ln(a,b,m);
    	ll ans=(b[n]*fac[n]%p+p)%p;
    	printf("%lld
    ",ans);
    	return 0;
    }
    
  • 相关阅读:
    css样式预处理器
    cookie&localStorage&sessionStorage
    《程序员面试金典》---读书笔记
    《Linux多线程服务端编程--使用muduo C++ 网络库》---读书笔记
    慢慢走上正轨
    padding 和 margin 不为人知的一面 ---(深度整理)
    html+css代码需要注意的地方(整理)
    前言
    MikTex 和 TexStudio 输入中文日文
    .htaccess 二级域名绑定子目录
  • 原文地址:https://www.cnblogs.com/ywwyww/p/8511071.html
Copyright © 2011-2022 走看看