码迷，mamicode.com
首页 > 其他好文 > 详细

有标号的DAG/强连通图计数

时间：2019-06-23 14:04:08      阅读：90      评论：0      收藏：0      [点我收藏+]

标签：php   inline   return   组合   for   ++   spl   不难   ons   

COGS索引

一堆神仙容斥+多项式……

有标号的DAG计数 I

考虑\(O(n^2)\)做法：设\(f_i\)表示总共有\(i\)个点的DAG数量，转移考虑枚举DAG上所有出度为\(0\)的点，剩下的点可以选择连向它，剩下的点之间也可以连边。

但是注意到这样子转移可能会存在剩下的点中有点没有出度的情况，考虑容斥解决：设枚举的出度为\(0\)的点的个数为\(i\)时的容斥系数为\(f_i\)，那么一个实际上存在\(x\)个出度为\(0\)的点的DAG的贡献就是\(\sum\limits_{i=1}^x \binom{x}{i} f_i = 1\)，不难由二项式定理知道\(f_i = (-1)^{i-1}\)

那么转移式就是\(f_i = \sum\limits_{j=1}^i \binom{i}{j} (-1)^{j-1} 2^{j(i-j)} f_{i-j}\)。

#include<bits/stdc++.h>
//this code is written by Itst
using namespace std;

const int MOD = 10007;
int dp[5003] , C[5003][5003] , poww2[6250003] , N;

int main(){
    freopen("DAG.in","r",stdin);
    freopen("DAG.out","w",stdout);
    cin >> N;
    poww2[0] = 1;
    for(int i = 1 ; i <= (N + 1) / 2 * ((N + 1) / 2) ; ++i)
        poww2[i] = (poww2[i - 1] << 1) % MOD;
    for(int i = 0 ; i <= N ; ++i){
        C[i][0] = 1;
        for(int j = 1 ; j <= i ; ++j)
            C[i][j] = (C[i - 1][j - 1] + C[i - 1][j]) % MOD;
    }
    dp[0] = 1;
    for(int i = 1 ; i <= N ; ++i)
        for(int j = 1 ; j <= i ; ++j)
            dp[i] = (dp[i] + ((j - 1) & 1 ? -1 : 1) * dp[i - j] * C[i][j] % MOD * poww2[j * (i - j)] % MOD + MOD) % MOD;
    cout << dp[N];
    return 0;
}

有标号的DAG计数 II

考虑使用多项式优化I中的做法。

一个前置芝士是使用组合数拆\(ij\)：\(ij = \binom{i}{2} + \binom{j+1}{2} - \binom{i-j}{2}\)

\[f_i = \sum\limits_{j=1}^i \binom{i}{j} (-1)^{j-1} 2^{j(i-j)} f_{i-j} \]

\[\frac{f_i}{i!} = \sum\limits_{j=1}^i (-1)^{j-1} 2^{\binom{i}{2} - \binom{j}{2} - \binom{i-j}{2}} \frac{f_{i-j}}{j!(i-j)!}\]

\[\frac{f_i}{i!2^\binom{i}{2}} = \sum\limits_{j=1}^i \frac{(-1)^{j-1}}{j!2^\binom{j}{2}} \frac{f_{i-j}}{(i-j)!2^\binom{i-j}{2}}\]

可以直接多项式求逆了。但是值得注意的一件事情是余项：因为\(j\)的下标从\(1\)开始，所以当\(i=0\)的时候，左式求出来为\(1\)，但是右式求出来为\(0\)。在多项式运算的时候记得补上这个余项。

#include<bits/stdc++.h>
//this code is written by Itst
using namespace std;

const int _ = (1 << 18) + 7 , MOD = 998244353;

int poww(long long a , long long b){
    int times = 1;
    while(b){
        if(b & 1) times = times * a % MOD;
        a = a * a % MOD; b >>= 1;
    }
    return times;
}

namespace poly{
    const int G = 3 , INV = 332748118;
    int dir[_] , need , invnd , A[_] , B[_];

    void init(int x){
        need = 1;
        while(need < x) need <<= 1;
        invnd = poww(need , MOD - 2);
        for(int i = 1 ; i < need ; ++i)
            dir[i] = (dir[i >> 1] >> 1) | (i & 1 ? need >> 1 : 0);
    }

    void NTT(int *arr , int tp){
        for(int i = 1 ; i < need ; ++i)
            if(i < dir[i])
                arr[i] ^= arr[dir[i]] ^= arr[i] ^= arr[dir[i]];
        for(int i = 1 ; i < need ; i <<= 1){
            int wn = poww(tp == 1 ? G : INV , MOD / i / 2);
            for(int j = 0 ; j < need ; j += i << 1){
                long long w = 1;
                for(int k = 0 ; k < i ; ++k , w = w * wn % MOD){
                    int x = arr[j + k] , y = arr[i + j + k] * w % MOD;
                    arr[j + k] = x + y >= MOD ? x + y - MOD : x + y;
                    arr[i + j + k] = x < y ? x + MOD - y : x - y;
                }
            }
        }
        if(tp != 1)
            for(int i = 0 ; i < need ; ++i)
                arr[i] = 1ll * arr[i] * invnd % MOD;
    }

#define clr(x) memset(x , 0 , sizeof(int) * need)
    void getInv(int *a , int *b , int len){
        if(len == 1){b[0] = poww(a[0] , MOD - 2); return;}
        getInv(a , b , ((len + 1) >> 1));
        memcpy(A , a , sizeof(int) * len); memcpy(B , b , sizeof(int) * len);
        init(len * 2 + 3); NTT(A , 1); NTT(B , 1);
        for(int i = 0 ; i < need ; ++i) A[i] = 1ll * A[i] * B[i] % MOD * B[i] % MOD;
        NTT(A , -1);
        for(int i = 0 ; i < len ; ++i)
            b[i] = (2ll * b[i] - A[i] + MOD) % MOD;
        clr(A); clr(B);
    }
}
using poly::getInv;
#define ch2(x) (1ll * x * (x - 1) / 2)
int F[_] , ans[_] , inv[_] , jc[_] , N;

void init(){
    jc[0] = 1;
    for(int i = 1 ; i <= N ; ++i) jc[i] = 1ll * jc[i - 1] * i % MOD;
    inv[N] = poww(jc[N] , MOD - 2);
    for(int i = N - 1 ; i >= 0 ; --i) inv[i] = inv[i + 1] * (i + 1ll) % MOD;
}

int main(){
    freopen("dag_count.in","r",stdin);
    freopen("dag_count.out","w",stdout);
    cin >> N; init();
    for(int i = 1 ; i <= N ; ++i) F[i] = (MOD + ((i - 1) & 1 ? 1ll : -1ll) * inv[i] * poww(poww(2 , ch2(i)) , MOD - 2) % MOD) % MOD;
    F[0] = 1; getInv(F , ans , N + 1);
    cout << 1ll * ans[N] * jc[N] % MOD * poww(2 , ch2(N)) % MOD;
    return 0;
}

有标号的DAG计数 III

在I和II中我们求出了可以不连通的DAG数量，而在这个部分我们强制要求DAG弱连通。

考虑用总的DAG数量减去不弱联通的DAG数量。设\(g_i\)表示点数为\(i\)的弱联通的DAG数量，总DAG数量就是I中求出的\(f_i\)，而对于不连通的DAG，它一定由若干个连通的DAG构成。那么我们考虑枚举\(1\)号点所在的弱联通DAG的大小，可以得到转移式：

\(g_i = f_i - \sum\limits_{j=1}^{i-1} f_jg_{i-j} \binom{i-1}{j-1}\)。直接转移复杂度\(O(n^2)\)。

#include<bits/stdc++.h>
//this code is written by Itst
using namespace std;

const int MOD = 10007;
int dp[5003] , g[5003] , C[5003][5003] , poww2[6250003] , N;

int main(){
    freopen("DAGIII.in","r",stdin);
    freopen("DAGIII.out","w",stdout);
    cin >> N;
    poww2[0] = 1;
    for(int i = 1 ; i <= (N + 1) / 2 * ((N + 1) / 2) ; ++i)
        poww2[i] = (poww2[i - 1] << 1) % MOD;
    for(int i = 0 ; i <= N ; ++i){
        C[i][0] = 1;
        for(int j = 1 ; j <= i ; ++j)
            C[i][j] = (C[i - 1][j - 1] + C[i - 1][j]) % MOD;
    }
    dp[0] = 1;
    for(int i = 1 ; i <= N ; ++i)
        for(int j = 1 ; j <= i ; ++j)
            dp[i] = (dp[i] + ((j - 1) & 1 ? -1 : 1) * dp[i - j] * C[i][j] % MOD * poww2[j * (i - j)] % MOD + MOD) % MOD;
    for(int i = 0 ; i <= N ; ++i){
        g[i] = dp[i];
        for(int j = 1 ; j < i ; ++j)
            g[i] = (g[i] - dp[j] * g[i - j] % MOD * C[i - 1][i - j - 1] % MOD + MOD) % MOD;
    }
    cout << g[N];
    return 0;
}

有标号的DAG计数 IV

有两种做法：

法一

考虑优化III中的递推式。

\[g_i = f_i - \sum\limits_{j=1}^{i-1} f_jg_{i-j} \binom{i-1}{j-1}\]

\[\frac{g_i}{(i-1)!} = \frac{f_i}{(i-1)!} - \sum\limits_{j=1}^{i-1} \frac{f_j}{(j-1)!}\frac{g_{i-j}}{(i-j)!}\]

多项式求逆即可。注意求和的下界和上界是\(1\)和\(i-1\)，所以\(g_0\)和\(f_0\)都要等于\(0\)才满足卷积式子。

#include<bits/stdc++.h>
//this code is written by Itst
using namespace std;

const int _ = (1 << 18) + 7 , MOD = 998244353;

int poww(long long a , long long b){
    int times = 1;
    while(b){
        if(b & 1) times = times * a % MOD;
        a = a * a % MOD; b >>= 1;
    }
    return times;
}

namespace poly{
    const int G = 3 , INV = 332748118;
    int dir[_] , need , invnd , A[_] , B[_];

    void init(int x){
        need = 1;
        while(need < x) need <<= 1;
        invnd = poww(need , MOD - 2);
        for(int i = 1 ; i < need ; ++i)
            dir[i] = (dir[i >> 1] >> 1) | (i & 1 ? need >> 1 : 0);
    }

    void NTT(int *arr , int tp){
        for(int i = 1 ; i < need ; ++i)
            if(i < dir[i])
                arr[i] ^= arr[dir[i]] ^= arr[i] ^= arr[dir[i]];
        for(int i = 1 ; i < need ; i <<= 1){
            int wn = poww(tp == 1 ? G : INV , MOD / i / 2);
            for(int j = 0 ; j < need ; j += i << 1){
                long long w = 1;
                for(int k = 0 ; k < i ; ++k , w = w * wn % MOD){
                    int x = arr[j + k] , y = arr[i + j + k] * w % MOD;
                    arr[j + k] = x + y >= MOD ? x + y - MOD : x + y;
                    arr[i + j + k] = x < y ? x + MOD - y : x - y;
                }
            }
        }
        if(tp != 1)
            for(int i = 0 ; i < need ; ++i)
                arr[i] = 1ll * arr[i] * invnd % MOD;
    }

#define clr(x) memset(x , 0 , sizeof(int) * need)
    void getInv(int *a , int *b , int len){
        if(len == 1){b[0] = poww(a[0] , MOD - 2); return;}
        getInv(a , b , ((len + 1) >> 1));
        memcpy(A , a , sizeof(int) * len); memcpy(B , b , sizeof(int) * len);
        init(len * 2 + 3); NTT(A , 1); NTT(B , 1);
        for(int i = 0 ; i < need ; ++i) A[i] = 1ll * A[i] * B[i] % MOD * B[i] % MOD;
        NTT(A , -1);
        for(int i = 0 ; i < len ; ++i)
            b[i] = (2ll * b[i] - A[i] + MOD) % MOD;
        clr(A); clr(B);
    }
}
using poly::getInv;
#define ch2(x) (1ll * x * (x - 1) / 2)
int F[_] , G[_] , ans[_] , tp[_] , inv[_] , jc[_] , N;

void init(){
    jc[0] = 1;
    for(int i = 1 ; i <= N ; ++i) jc[i] = 1ll * jc[i - 1] * i % MOD;
    inv[N] = poww(jc[N] , MOD - 2);
    for(int i = N - 1 ; i >= 0 ; --i) inv[i] = inv[i + 1] * (i + 1ll) % MOD;
}

int main(){
    freopen("dagIV.in","r",stdin);
    freopen("dagIV.out","w",stdout);
    cin >> N; init();
    for(int i = 1 ; i <= N ; ++i) F[i] = (MOD + ((i - 1) & 1 ? 1ll : -1ll) * inv[i] * poww(poww(2 , ch2(i)) , MOD - 2) % MOD) % MOD;
    F[0] = 1; getInv(F , G , N + 1);
    for(int i = 1 ; i <= N ; ++i)
        G[i] = 1ll * G[i] * jc[i] % MOD * poww(2 , ch2(i)) % MOD;
    for(int i = 0 ; i <= N ; ++i)
        tp[i] = 1ll * G[i] * inv[i] % MOD;
    getInv(tp , ans , N + 1); poly::init(2 * N + 2);
    for(int i = 1 ; i <= N ; ++i)
        G[i] = 1ll * G[i] * inv[i - 1] % MOD;
    G[0] = 0;
    poly::NTT(G , 1); poly::NTT(ans , 1);
    for(int i = 0 ; i < poly::need ; ++i)
        ans[i] = 1ll * ans[i] * G[i] % MOD;
    poly::NTT(ans , 0);
    cout << 1ll * ans[N] * jc[N - 1] % MOD;
    return 0;
}

法二

注意到III中的一句话：一个不一定连通的DAG一定由若干个连通的DAG构成，这等价于一个不一定连通的DAG是若干个连通的DAG的有标号集合。计算有标号集合可以考虑多项式Exp，那么如果设连通的DAG的数量的指数型生成函数为\(F\)，不一定连通的DAG的数量的指数型生成函数为\(G\)，那么\(G = e^F\)，即\(F = ln\ G\)。多项式求Ln即可。

#include<bits/stdc++.h>
//this code is written by Itst
using namespace std;

const int _ = (1 << 18) + 7 , MOD = 998244353;

int poww(long long a , long long b){
    int times = 1;
    while(b){
        if(b & 1) times = times * a % MOD;
        a = a * a % MOD; b >>= 1;
    }
    return times;
}

namespace poly{
    const int G = 3 , INV = 332748118;
    int dir[_] , need , invnd , A[_] , B[_] , C[_];

    void init(int x){
        need = 1;
        while(need < x) need <<= 1;
        invnd = poww(need , MOD - 2);
        for(int i = 1 ; i < need ; ++i)
            dir[i] = (dir[i >> 1] >> 1) | (i & 1 ? need >> 1 : 0);
    }

    void NTT(int *arr , int tp){
        for(int i = 1 ; i < need ; ++i)
            if(i < dir[i])
                arr[i] ^= arr[dir[i]] ^= arr[i] ^= arr[dir[i]];
        for(int i = 1 ; i < need ; i <<= 1){
            int wn = poww(tp == 1 ? G : INV , MOD / i / 2);
            for(int j = 0 ; j < need ; j += i << 1){
                long long w = 1;
                for(int k = 0 ; k < i ; ++k , w = w * wn % MOD){
                    int x = arr[j + k] , y = arr[i + j + k] * w % MOD;
                    arr[j + k] = x + y >= MOD ? x + y - MOD : x + y;
                    arr[i + j + k] = x < y ? x + MOD - y : x - y;
                }
            }
        }
        if(tp != 1)
            for(int i = 0 ; i < need ; ++i)
                arr[i] = 1ll * arr[i] * invnd % MOD;
    }

#define clr(x) memset(x , 0 , sizeof(int) * need)
    void getInv(int *a , int *b , int len){
        if(len == 1){b[0] = poww(a[0] , MOD - 2); return;}
        getInv(a , b , ((len + 1) >> 1));
        memcpy(A , a , sizeof(int) * len); memcpy(B , b , sizeof(int) * len);
        init(len * 2 + 3); NTT(A , 1); NTT(B , 1);
        for(int i = 0 ; i < need ; ++i) A[i] = 1ll * A[i] * B[i] % MOD * B[i] % MOD;
        NTT(A , -1);
        for(int i = 0 ; i < len ; ++i)
            b[i] = (2ll * b[i] - A[i] + MOD) % MOD;
        clr(A); clr(B);
    }

    void getDis(int *a , int *b , int len){
        for(int i = 0 ; i < len - 1 ; ++i)
            b[i] = a[i + 1] * (i + 1ll) % MOD;
        b[len - 1] = 0;
    }

    void getInt(int *a , int *b , int len){
        for(int i = 1 ; i <= len ; ++i)
            b[i] = 1ll * a[i - 1] * poww(i , MOD - 2) % MOD;
    }

    void getLn(int *a , int *b , int len){
        getInv(a , C , len); getDis(a , A , len);
        init(2 * len + 1); NTT(A , 1); NTT(C , 1);
        for(int i = 0 ; i < need ; ++i)
            A[i] = 1ll * A[i] * C[i] % MOD;
        NTT(A , -1);
        getInt(A , b , len - 1); clr(A); clr(C);
    }
}
using poly::getInv; using poly::getLn;
#define ch2(x) (1ll * x * (x - 1) / 2)
int F[_] , G[_] , ans[_] , tp[_] , inv[_] , jc[_] , N;

void init(){
    jc[0] = 1;
    for(int i = 1 ; i <= N ; ++i) jc[i] = 1ll * jc[i - 1] * i % MOD;
    inv[N] = poww(jc[N] , MOD - 2);
    for(int i = N - 1 ; i >= 0 ; --i) inv[i] = inv[i + 1] * (i + 1ll) % MOD;
}

int main(){
    freopen("dagIV.in","r",stdin);
    freopen("dagIV.out","w",stdout);
    cin >> N; init();
    for(int i = 1 ; i <= N ; ++i) F[i] = (MOD + ((i - 1) & 1 ? 1ll : -1ll) * inv[i] * poww(poww(2 , ch2(i)) , MOD - 2) % MOD) % MOD;
    F[0] = 1; getInv(F , G , N + 1);
    for(int i = 1 ; i <= N ; ++i)
        G[i] = 1ll * G[i] * poww(2 , ch2(i)) % MOD;
    getLn(G , ans , N + 1);
    cout << 1ll * ans[N] * jc[N] % MOD;
    return 0;
}

有标号的强连通图计数 I

感觉比上面的难不少但是评级却更低是什么鬼

一个任意的有向图在缩点时候都可以得到一个DAG，那么我们可以仍然考虑DAG计数I中的做法，即枚举缩点之后的图中出度为\(0\)的点由哪些点构成。但是值得注意的一件事情是在DAG计数I中，容斥系数与出度为\(0\)的点数相关，那么在强连通图计数中，这个容斥系数应该与选择的点在缩点之后构成的强连通分量的个数相关，相比于DAG计数来说这个是最为棘手的。

考虑设\(f_i\)表示\(i\)个点的强连通图数量，$g_i = \sum\limits_{j=1}^i (-1)^{j-1} \times $\(i\)个点构成\(j\)个强连通分量的方案数，\(h_i = 2^{i(i-1)}\)表示\(i\)个点的有向图数量。不难发现\(g_i\)就是把容斥系数和方案数放在了一起。

那么枚举缩点之后出度为\(0\)的强连通分量由哪些点构成，可以得到转移式：\(f_i = h_i - \sum\limits_{j=1}^i \binom{i}{j} h_{i-j} g_{j} 2^{j(i-j)} + f_i\)。值得注意的是最后加上的\(f_i\)，这是因为当\(j=i\)时，\(g_i\)中包含了\(i\)个点构成同一个强连通分量的方案数，所以要把它补回来。

由上面的式子变形一下就可以得到\(g_i\)的转移式：\(g_i = h_i - \sum\limits_{j=1}^{i-1} \binom{i}{j} h_{i-j} g_{j} 2^{j(i-j)}\)。

然后考虑\(f\)和\(g\)的关系。按照上面的定义，\(g_i\)表示的是\(i\)个点构成强连通分量的方案总和，其中奇数个强连通分量的方案贡献为\(1\)，偶数个强连通分量的贡献为\(-1\)。那么我们考虑枚举\(1\)号点所在的强连通分量，可以得到转移式：\(g_i = f_i - \sum\limits_{j=1}^{i-1} f_j g_{i-j} \binom{i-1}{j-1}\)，即\(f_i = g_i + \sum\limits_{j=1}^{i-1} f_j g_{i-j} \binom{i-1}{j-1}\)。

暴力求出\(f,g\)，复杂度\(O(n^2)\)。

#include<bits/stdc++.h>
//this code is written by Itst
using namespace std;

const int _ = 1003 , MOD = 10007;

int poww(long long a , long long b){
    int times = 1;
    while(b){
        if(b & 1) times = times * a % MOD;
        a = a * a % MOD; b >>= 1;
    }
    return times;
}

int F[_] , G[_] , inv[_] , jc[_] , poww2[_ * _] , N;

void init(){
    jc[0] = 1;
    for(int i = 1 ; i <= N ; ++i) jc[i] = jc[i - 1] * i % MOD;
    inv[N] = poww(jc[N] , MOD - 2);
    for(int i = N - 1 ; i >= 0 ; --i) inv[i] = inv[i + 1] * (i + 1) % MOD;
    poww2[0] = 1;
    for(int i = 1 ; i <= N * N ; ++i)
        poww2[i] = poww2[i - 1] * 2 % MOD;
}

int binom(int a , int b){return a < b ? 0 : jc[a] * inv[b] % MOD * inv[a - b] % MOD;}

int main(){
    freopen("QAQ_strong_one.in","r",stdin);
    freopen("QAQ_strong_one.out","w",stdout);
    cin >> N; init();
    for(int i = 1 ; i <= N ; ++i){
        G[i] = poww2[i * (i - 1)];
        for(int j = 1 ; j < i ; ++j)
            G[i] = (G[i] - binom(i , j) * poww2[j * (i - j)] % MOD * poww2[(i - j) * (i - j - 1)] % MOD * G[j] % MOD + MOD) % MOD;
    }
    for(int i = 1 ; i <= N ; ++i){
        F[i] = G[i];
        for(int j = 1 ; j < i ; ++j)
            F[i] = (F[i] + F[j] * G[i - j] % MOD * binom(i - 1 , j - 1)) % MOD;
    }
    cout << F[N];
    return 0;
}

有标号的强连通图计数 II

第一部分：多项式优化求\(g\)。和DAG计数II差不多，推一下式子可以得到\(\frac{g_n}{n!2^\binom{n}{2}} = \frac{h_n}{n!2^\binom{n}{2}} - \sum\limits_{j=1}^{n-1} \frac{h_{n-j}}{(n-j)! 2^\binom{n-j}{2}} \frac{g_j}{j!2^\binom{j}{2}}\)，仍然是多项式求逆，仍然需要注意边界。

第二部分：多项式优化求\(f\)。基本思路和DAG计数IV是一样的，有拆式子然后多项式求逆的方法，也可以考虑性质然后多项式Ln解决。

多项式Ln代码

#include<bits/stdc++.h>
//this code is written by Itst
using namespace std;

const int _ = (1 << 18) + 7 , MOD = 998244353;

int poww(long long a , long long b){
    int times = 1;
    while(b){
        if(b & 1) times = times * a % MOD;
        a = a * a % MOD; b >>= 1;
    }
    return times;
}

namespace poly{
    const int G = 3 , INV = 332748118;
    int dir[_] , need , invnd , A[_] , B[_] , C[_];

    void init(int x){
        need = 1;
        while(need < x) need <<= 1;
        invnd = poww(need , MOD - 2);
        for(int i = 1 ; i < need ; ++i)
            dir[i] = (dir[i >> 1] >> 1) | (i & 1 ? need >> 1 : 0);
    }

    void NTT(int *arr , int tp){
        for(int i = 1 ; i < need ; ++i)
            if(i < dir[i])
                arr[i] ^= arr[dir[i]] ^= arr[i] ^= arr[dir[i]];
        for(int i = 1 ; i < need ; i <<= 1){
            int wn = poww(tp == 1 ? G : INV , MOD / i / 2);
            for(int j = 0 ; j < need ; j += i << 1){
                long long w = 1;
                for(int k = 0 ; k < i ; ++k , w = w * wn % MOD){
                    int x = arr[j + k] , y = arr[i + j + k] * w % MOD;
                    arr[j + k] = x + y >= MOD ? x + y - MOD : x + y;
                    arr[i + j + k] = x < y ? x + MOD - y : x - y;
                }
            }
        }
        if(tp != 1)
            for(int i = 0 ; i < need ; ++i)
                arr[i] = 1ll * arr[i] * invnd % MOD;
    }

#define clr(x) memset(x , 0 , sizeof(int) * need)
    void getInv(int *a , int *b , int len){
        if(len == 1){b[0] = poww(a[0] , MOD - 2); return;}
        getInv(a , b , ((len + 1) >> 1));
        memcpy(A , a , sizeof(int) * len); memcpy(B , b , sizeof(int) * len);
        init(len * 2 + 3); NTT(A , 1); NTT(B , 1);
        for(int i = 0 ; i < need ; ++i) A[i] = 1ll * A[i] * B[i] % MOD * B[i] % MOD;
        NTT(A , -1);
        for(int i = 0 ; i < len ; ++i)
            b[i] = (2ll * b[i] - A[i] + MOD) % MOD;
        clr(A); clr(B);
    }

    void getDis(int *a , int *b , int len){
        for(int i = 0 ; i < len - 1 ; ++i)
            b[i] = a[i + 1] * (i + 1ll) % MOD;
        b[len - 1] = 0;
    }

    void getInt(int *a , int *b , int len){
        for(int i = 1 ; i <= len ; ++i)
            b[i] = 1ll * a[i - 1] * poww(i , MOD - 2) % MOD;
    }

    void getLn(int *a , int *b , int len){
        getInv(a , C , len); getDis(a , A , len);
        init(2 * len + 1); NTT(A , 1); NTT(C , 1);
        for(int i = 0 ; i < need ; ++i)
            A[i] = 1ll * A[i] * C[i] % MOD;
        NTT(A , -1);
        getInt(A , b , len - 1); clr(A); clr(C);
    }
}
using poly::getInv; using poly::getLn;
#define ch2(x) (1ll * x * (x - 1) / 2)
int F[_] , W[_] , H[_] , inv[_] , jc[_] , N;

void init(){
    jc[0] = 1;
    for(int i = 1 ; i <= N ; ++i) jc[i] = 1ll * jc[i - 1] * i % MOD;
    inv[N] = poww(jc[N] , MOD - 2);
    for(int i = N - 1 ; i >= 0 ; --i) inv[i] = inv[i + 1] * (i + 1ll) % MOD;
}

int main(){
    freopen("QAQ_strongly_two.in","r",stdin);
    freopen("QAQ_strongly_two.out","w",stdout);
    cin >> N; init();
    for(int i = 1 ; i <= N ; ++i)
        H[i] = 1ll * poww(2 , ch2(i)) * inv[i] % MOD;
    H[0] = 1; getInv(H , W , N + 1); H[0] = 0;
    poly::init(2 * N + 2); poly::NTT(H , 1); poly::NTT(W , 1);
    for(int i = 0 ; i < poly::need ; ++i)
        H[i] = 1ll * H[i] * W[i] % MOD;
    poly::NTT(H , -1); H[0] = 1;
    for(int i = 1 ; i <= N ; ++i)
        H[i] = MOD - 1ll * H[i] * poww(2 , ch2(i)) % MOD;
    getLn(H , F , N + 1);
    cout << MOD - 1ll * F[N] * jc[N] % MOD;
    return 0;
}

多项式求逆代码

#include<bits/stdc++.h>
//this code is written by Itst
using namespace std;

const int _ = (1 << 18) + 7 , MOD = 998244353;

int poww(long long a , long long b){
    int times = 1;
    while(b){
        if(b & 1) times = times * a % MOD;
        a = a * a % MOD; b >>= 1;
    }
    return times;
}

namespace poly{
    const int G = 3 , INV = 332748118;
    int dir[_] , need , invnd , A[_] , B[_] , C[_];

    void init(int x){
        need = 1;
        while(need < x) need <<= 1;
        invnd = poww(need , MOD - 2);
        for(int i = 1 ; i < need ; ++i)
            dir[i] = (dir[i >> 1] >> 1) | (i & 1 ? need >> 1 : 0);
    }

    void NTT(int *arr , int tp){
        for(int i = 1 ; i < need ; ++i)
            if(i < dir[i])
                arr[i] ^= arr[dir[i]] ^= arr[i] ^= arr[dir[i]];
        for(int i = 1 ; i < need ; i <<= 1){
            int wn = poww(tp == 1 ? G : INV , MOD / i / 2);
            for(int j = 0 ; j < need ; j += i << 1){
                long long w = 1;
                for(int k = 0 ; k < i ; ++k , w = w * wn % MOD){
                    int x = arr[j + k] , y = arr[i + j + k] * w % MOD;
                    arr[j + k] = x + y >= MOD ? x + y - MOD : x + y;
                    arr[i + j + k] = x < y ? x + MOD - y : x - y;
                }
            }
        }
        if(tp != 1)
            for(int i = 0 ; i < need ; ++i)
                arr[i] = 1ll * arr[i] * invnd % MOD;
    }

#define clr(x) memset(x , 0 , sizeof(int) * need)
    void getInv(int *a , int *b , int len){
        if(len == 1){b[0] = poww(a[0] , MOD - 2); return;}
        getInv(a , b , ((len + 1) >> 1));
        memcpy(A , a , sizeof(int) * len); memcpy(B , b , sizeof(int) * len);
        init(len * 2 + 3); NTT(A , 1); NTT(B , 1);
        for(int i = 0 ; i < need ; ++i) A[i] = 1ll * A[i] * B[i] % MOD * B[i] % MOD;
        NTT(A , -1);
        for(int i = 0 ; i < len ; ++i)
            b[i] = (2ll * b[i] - A[i] + MOD) % MOD;
        clr(A); clr(B);
    }
}
using poly::getInv;
#define ch2(x) (1ll * x * (x - 1) / 2)
int F[_] , W[_] , H[_] , inv[_] , jc[_] , N;

void init(){
    jc[0] = 1;
    for(int i = 1 ; i <= N ; ++i) jc[i] = 1ll * jc[i - 1] * i % MOD;
    inv[N] = poww(jc[N] , MOD - 2);
    for(int i = N - 1 ; i >= 0 ; --i) inv[i] = inv[i + 1] * (i + 1ll) % MOD;
}

int main(){
    freopen("QAQ_strongly_two.in","r",stdin);
    freopen("QAQ_strongly_two.out","w",stdout);
    cin >> N; init();
    for(int i = 1 ; i <= N ; ++i)
        H[i] = 1ll * poww(2 , ch2(i)) * inv[i] % MOD;
    H[0] = 1; getInv(H , W , N + 1); H[0] = 0;
    poly::init(2 * N + 2); poly::NTT(H , 1); poly::NTT(W , 1);
    for(int i = 0 ; i < poly::need ; ++i)
        H[i] = 1ll * H[i] * W[i] % MOD;
    poly::NTT(H , -1); H[0] = 1;
    for(int i = 1 ; i <= N ; ++i)
        H[i] = MOD - 1ll * H[i] * poww(2 , ch2(i)) % MOD;
    getInv(H , F , N + 1); H[0] = 0;
    for(int i = 1 ; i <= N ; ++i)
        H[i] = 1ll * (MOD - H[i]) * i % MOD;
    poly::init(2 * N + 2); poly::NTT(H , 1); poly::NTT(F , 1);
    for(int i = 0 ; i < poly::need ; ++i)
        H[i] = 1ll * H[i] * F[i] % MOD;
    poly::NTT(H , -1);
    cout << 1ll * H[N] * jc[N - 1] % MOD;
    return 0;
}

有标号的DAG/强连通图计数

标签：php   inline   return   组合   for   ++   spl   不难   ons   

原文地址：https://www.cnblogs.com/Itst/p/11072400.html
踩
(0)
赞
(0)
   
举报
评论 一句话评论（0）
登录后才能评论！
分享档案
更多>
2021年07月29日 (22)
2021年07月28日 (40)
2021年07月27日 (32)
2021年07月26日 (79)
2021年07月23日 (29)
2021年07月22日 (30)
2021年07月21日 (42)
2021年07月20日 (16)
2021年07月19日 (90)
2021年07月16日 (35)
周排行
mamicode.com排行更多图片
更多
友情链接
兰亭集智   国之画   百度统计   站长统计   阿里云   chrome插件   新版天听网
关于我们 - 联系我们 - 留言反馈
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码！