zoukankan      html  css  js  c++  java
  • 【bzoj3456】城市规划 容斥原理+NTT+多项式求逆

    题目描述

    求出n个点的简单(无重边无自环)无向连通图数目mod 1004535809(479 * 2 ^ 21 + 1).

    输入

    仅一行一个整数n(<=130000)

    输出

    仅一行一个整数, 为方案数 mod 1004535809.

    样例输入

    3

    样例输出

    4


    题解

    容斥原理+NTT+多项式求逆

    设 $f_i$ 表示 $i$ 个点的简单无向连通图的数目,$g_i$ 表示 $i$ 个点的简单无向图的数目。

    根据定义得 $g_i=2^{frac{n(n-1}2}$ 。

    对于 $f_i$ ,考虑容斥,用 $g_i$ 减去不连通的方案数。枚举不连通图中1号点所在连通块大小 $j$ ,则有:

    $$f_i=g_i-sumlimits_{j=1}^{i-1}C_{i-1}^{j-1}f_jg_{i-j}$$ 

    将组合数展开,得:

    $$f_i=g_i-sumlimits_{j=1}^{i-1}frac{(i-1)!}{(j-1)!(i-j)!}f_jg_{i-j}$$ 

    两边同时除以 $(i-1)!$ ,整理得:

    $$frac{f_i}{(i-1)!}=sumlimits_{j=1}^{i-1}frac{f_j}{(j-1)!}frac{g_{i-j}}{(i-j)!}$$ 

    设:

    $$F(x)=sumlimits_{i=1}^{infty}frac{f_i}{(i-1)!} \ G(x)=sumlimits_{i=1}^{infty}frac{g_i}{(i-1)!} \ H(x)=sumlimits_{i=1}^{infty}frac{g_i}{i!}$$ 

    注意到 $F(x)$ 、$H(x)$ 的常数项均为0,因此后面的那个式子相当于 $sumlimits_{j=0}^iF(x)[j]H(x)[i-j]$ ,是一个卷积的形式。

    因此有:

    $$F(x)=G(x)-F(x) imes H(x)$$ 

    化简得:

    $$F(x)=frac{G(x)}{H(x)+1}$$ 

    剩下的就好办了,根据定义求出 $G(x)$ 和 $H(x)+1$ ,使用多项式求逆求出 $H(x)+1$ 的逆,再与 $G(x)$ 求乘法即可得到 $F(x)$ 。

    最后的答案就是 $F(x)[n] imes(n-1)!$ 。

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

    #include <cstdio>
    #include <cstring>
    #include <algorithm>
    #define N 262200
    #define mod 1004535809
    using namespace std;
    typedef long long ll;
    ll A[N] , B[N] , C[N] , t[N] , fac[N >> 1];
    inline ll pow(ll x , ll y)
    {
        ll ans = 1;
        while(y)
        {
            if(y & 1) ans = ans * x % mod;
            x = x * x % mod , y >>= 1;
        }
        return ans;
    }
    void ntt(ll *a , int n , int flag)
    {
        int i , j , k;
        for(i = k = 0 ; i < n ; i ++ )
        {
            if(i > k) swap(a[i] , a[k]);
            for(j = n >> 1 ; (k ^= j) < j ; j >>= 1);
        }
        for(k = 2 ; k <= n ; k <<= 1)
        {
            ll wn = pow(3 , (mod - 1) / k);
            if(flag == -1) wn = pow(wn , mod - 2);
            for(i = 0 ; i < n ; i += k)
            {
                ll w = 1 , t;
                for(j = i ; j < i + (k >> 1) ; j ++ , w = w * wn % mod)
                    t = w * a[j + (k >> 1)] % mod , a[j + (k >> 1)] = (a[j] - t + mod) % mod , a[j] = (a[j] + t) % mod;
            }
        }
        if(flag == -1)
        {
            k = pow(n , mod - 2);
            for(i = 0 ; i < n ; i ++ ) a[i] = a[i] * k % mod;
        }
    }
    void inv(ll *a , ll *b , int n)
    {
        if(n == 1)
        {
            b[0] = 1;
            return;
        }
        int i;
        inv(a , b , n >> 1);
        memcpy(t , a , sizeof(ll) * n);
        ntt(t , n << 1 , 1);
        ntt(b , n << 1 , 1);
        for(i = 0 ; i < n << 1 ; i ++ ) b[i] = b[i] * (2 - t[i] * b[i] % mod + mod) % mod;
        ntt(b , n << 1 , -1);
        memset(b + n , 0 , sizeof(ll) * n);
    }
    int main()
    {
        int n , i , len = 1;
        scanf("%d" , &n);
        fac[0] = A[0] = 1;
        for(i = 1 ; i <= n ; i ++ )
        {
            fac[i] = fac[i - 1] * i % mod;
            A[i] = pow(2 , (ll)i * (i - 1) / 2) * pow(fac[i] , mod - 2) % mod;
            B[i] = pow(2 , (ll)i * (i - 1) / 2) * pow(fac[i - 1] , mod - 2) % mod;
        }
        while(len <= n) len <<= 1;
        inv(A , C , len);
        ntt(B , len , 1);
        ntt(C , len , 1);
        for(i = 0 ; i < len ; i ++ ) B[i] = B[i] * C[i] % mod;
        ntt(B , len , -1);
        printf("%lld
    " , B[n] * fac[n - 1] % mod);
        return 0;
    }
    

     

  • 相关阅读:
    如何用Python实现网络请求库中的UR解析器,面试必学
    为什么有人说 Python 多线程是鸡肋?
    router-view 与 动态组件 区别
    keep-alive
    vue 客户端渲染和服务端渲染
    js 数组对象深拷贝
    vue template标签
    Jquery中的日历插件
    HTML5中的canvas基本概念及绘图
    HTML5中的音视频处理
  • 原文地址:https://www.cnblogs.com/GXZlegend/p/8596008.html
Copyright © 2011-2022 走看看