HDU 5552 (CDQ分治 NTT)

最新推荐文章于 2024-01-26 02:49:59 发布

morejarphone

最新推荐文章于 2024-01-26 02:49:59 发布

阅读量1.3k

点赞数

分类专栏： FFT && NTT 分治

本文链接：https://blog.csdn.net/morejarphone/article/details/52432921

版权

分治同时被 2 个专栏收录

12 篇文章 0 订阅

订阅专栏

FFT && NTT

9 篇文章 0 订阅

订阅专栏

题意：n个点有环连通图计数。每条边可以染成m种颜色。

定义一些数组：
f[n]: n个点的连通图；
g[n]:n个点的图；
h[n]:n个点的树。
后两者计算很简单，每条边都有 $m+1$ 种选择（不取这条边或者m种颜色），一共有 $\frac{n*(n-1)}{2}$ 条边，所以 $g[n]=(m+1)^\frac{n*(n-1)}{2}$ ；n个点的树根据prufer数列一个树对应一个n-2长度的序列，而一个树有n-1条边，故 $h[n]=n^{n-2}*m^{n-1}$ .

对于f数组， $联通图的个数=图的总数-非联通图的个数$ ，所以枚举1这个点所在的联通块的大小：

f [n] = g [n] - \sum i = 1 n - 1 C i - 1 n - 1 f [i] g [n - i] = g [n] - \sum i = 1 n - 1 ( n - 1 ) ! ( i - 1 ) ! * ( n - i ) ! = g [n] - (n - 1)! \sum i = 1 n - 1 f [ i ] ( i - 1 ) ! * g [ n - i ] ( n - i ) !

$\begin{align} f[n]&=g[n]-\sum_{i=1}^{n-1}C_{n-1}^{i-1}f[i]g[n-i] \\ &= g[n]-\sum_{i=1}^{n-1}\frac{(n-1)!}{(i-1)!*(n-i)!}\\ &=g[n]-(n-1)!\sum_{i=1}^{n-1}\frac{f[i]}{(i-1)!} * \frac{g[n-i]}{(n-i)!}\\ \end{align}$
这个东西是一个卷积的形式，用fft或者ntt加速就好了。

那么最后的答案就是 $n个点的连通图-n个点的树$ ，也就是 $ans=f[n]-h[n]$ 。

然后fft开心的跪了精度，所以用ntt来避免精度误差，原根的算出来是106。

FFT：（精度会炸）

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <iostream>
#include <cmath>
#include <map>
#include <vector>
#include <stack>
using namespace std;
#define mod 152076289
#define pi acos(-1.0)
#define maxn 40005

struct plex {
    long double x, y;
    plex (long double _x = 0.0, long double _y = 0.0) : x (_x), y (_y) {}
    plex operator + (const plex &a) const {
        return plex (x+a.x, y+a.y);
    }
    plex operator - (const plex &a) const {
        return plex (x-a.x, y-a.y);
    }
    plex operator * (const plex &a) const {
        return plex (x*a.x-y*a.y, x*a.y+y*a.x);
    }
}x1[maxn], x2[maxn];

void change (plex *y, int len) {
    int i, j, k;
    for(i = 1, j = len / 2; i < len - 1; i++) {
        if (i < j) swap(y[i], y[j]);
        k = len / 2;
        while (j >= k) {
            j -= k;
            k /= 2;
        }
        if (j < k) j += k;
    }
}

void fft(plex y[],int len,int on)
{
    change(y,len);
    for(int h = 2; h <= len; h <<= 1)
    {
        plex wn(cos(-on*2*pi/h),sin(-on*2*pi/h));
        for(int j = 0;j < len;j+=h)
        {
            plex w(1,0);
            for(int k = j;k < j+h/2;k++)
            {
                plex u = y[k];
                plex t = w*y[k+h/2];
                y[k] = u+t;
                y[k+h/2] = u-t;
                w = w*wn;
            }
        }
    }
    if(on == -1)
        for(int i = 0;i < len;i++)
            y[i].x /= len;
}

long long qpow (long long a, long long b) {  
    long long ret=1;  
    while (b) {  
        if (b&1) ret = (ret*a)%mod;  
        a = (a*a)%mod;  
        b >>= 1;  
    }  
    return ret;  
}  

long long n, m;
long long f[maxn], g[maxn], h[maxn];
//n个点的连通图 n个点的图 n个点的树
long long c[maxn], fac[maxn], rev[maxn];

void solve (int l, int r) {
    if (l == r) {
        f[l] += g[l];
        f[l] %= mod;
        return;
    }
    int mid = (l+r) >> 1;
    solve (l, mid);
    int len = 1;
    while (len <= r-l+1) {len <<= 1;}
    for (int i = 0; i < len; i++) {
        x1[i] = x2[i] = plex (0, 0);
    }
    for (int i = l; i <= mid; i++) {
        long long tmp = f[i]*rev[i-1]%mod;
        x1[i-l] = plex (tmp, 0);
    }
    for (int i = 1; i <= r-l; i++) {
        long long tmp = g[i]*rev[i]%mod;
        x2[i-1] = plex (tmp, 0);
    }
    fft (x1, len, 1), fft (x2, len, 1);
    for (int i = 0; i < len; i++) x1[i] = x1[i] * x2[i];
    fft (x1, len, -1);
    for (int i = mid+1; i <= r; i++) {
        f[i] -= (long long) (x1[i-l-1].x + 0.5) %mod *fac[i-1] %mod;
        (f[i] += mod) %= mod;
    }
    //cout << l << " " << r << ": ";
    //for (int i = 1; i <= n; i++) cout << f[i] << " "; cout << endl;
    solve (mid+1, r);
}

int main () {
    int t, kase = 0;
    fac[0] = 1;
    for (int i = 1; i < maxn; i++) fac[i] = fac[i-1] * i % mod;
    rev[maxn-1] = qpow (fac[maxn-1], mod - 2);
    for (int i = maxn-2; i >= 0; i--) rev[i] = rev[i+1] * (i+1) % mod;
    scanf ("%d", &t);
    while (t--) {
        scanf ("%lld%lld", &n, &m);
        memset (f, 0, sizeof f);
        c[0] = 1;
        for (int i = 1; i < n; i++) {
            c[i] = c[i-1]*(n-1-i+1)%mod*rev[i]%mod;
        }
        for (int i = 1; i <= n; i++) {
            g[i] = qpow (1+m, i*(i-1)/2);
        }
        //for (int i = 1; i <= n; i++) cout << c[i] << " "; cout << endl;
        //for (int i = 1; i <= n; i++) cout << g[i] << " "; cout << endl;
        solve (1, n);
        //for (int i = 1; i <= n; i++) cout << f[i] << " "; cout << endl;
        long long ans = f[n]; 
        long long tmp = qpow (n, n-2) * qpow (m, n-1) % mod;
        ans = (ans - tmp + mod) % mod;
        printf ("Case #%d: %lld\n", ++kase, ans);
    }
    return 0;
}

NTT：

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <iostream>
#include <cmath>
#include <map>
#include <vector>
#include <stack>
using namespace std;
#define mod 152076289
#define G 106
#define maxn 40005

long long x1[maxn], x2[maxn];
long long qpow (long long a, long long b) {  
    long long ret=1;  
    while (b) {  
        if (b&1) ret = (ret*a)%mod;  
        a = (a*a)%mod;  
        b >>= 1;  
    }  
    return ret;  
}  

void change (long long *y, int len) {
    int i, j, k;
    for(i = 1, j = len / 2; i < len - 1; i++) {
        if (i < j) swap(y[i], y[j]);
        k = len / 2;
        while (j >= k) {
            j -= k;
            k /= 2;
        }
        if (j < k) j += k;
    }
}

void ntt (long long *y, int len, int on) {
    change (y, len);
    int id = 0;
    for(int h = 2; h <= len; h <<= 1) {
        id++;
        long long wn = qpow (G, (mod - 1) / (1<<id));
        for(int j = 0; j < len; j += h) {
            long long w = 1;
            for(int k = j; k < j + h / 2; k++) {
                long long u = y[k] % mod;
                long long t = w * (y[k + h / 2] % mod) % mod;
                y[k] = (u + t) % mod;
                y[k + h / 2] = ((u - t) % mod + mod) % mod;
                w = w * wn % mod;
            }
        }
    }
    if (on == -1) {
        for (int i = 1; i < len / 2; i++)
            swap (y[i], y[len - i]);
        long long inv = qpow(len, mod - 2);
        for(int i = 0; i < len; i++)
            y[i] = y[i] % mod * inv % mod;
    }
}

long long n, m;
long long f[maxn], g[maxn], h[maxn];
//n个点的连通图 n个点的图 n个点的树
long long c[maxn], fac[maxn], rev[maxn];

void solve (int l, int r) {
    if (l == r) {
        f[l] += g[l];
        f[l] %= mod;
        return;
    }
    int mid = (l+r) >> 1;
    solve (l, mid);
    int len = 1;
    while (len <= r-l+1) {len <<= 1;}
    for (int i = 0; i < len; i++) {
        x1[i] = x2[i] = 0;
    }
    for (int i = l; i <= mid; i++) {
        x1[i-l] = f[i]*rev[i-1]%mod;
    }
    for (int i = 1; i <= r-l; i++) {
        x2[i-1] = g[i]*rev[i]%mod;
    }
    ntt (x1, len, 1), ntt (x2, len, 1);
    for (int i = 0; i < len; i++) x1[i] = x1[i] * x2[i] % mod;
    ntt (x1, len, -1);
    for (int i = mid+1; i <= r; i++) {
        f[i] -= x1[i-l-1] %mod *fac[i-1] %mod;
        (f[i] += mod) %= mod;
    }
    solve (mid+1, r);
}

int main () {
    int t, kase = 0;
    fac[0] = 1;
    for (int i = 1; i < maxn; i++) fac[i] = fac[i-1] * i % mod;
    rev[maxn-1] = qpow (fac[maxn-1], mod - 2);
    for (int i = maxn-2; i >= 0; i--) rev[i] = rev[i+1] * (i+1) % mod;
    scanf ("%d", &t);
    while (t--) {
        scanf ("%lld%lld", &n, &m);
        memset (f, 0, sizeof f);
        for (int i = 1; i <= n; i++) {
            g[i] = qpow (1+m, i*(i-1)/2);
        }
        solve (1, n);
        long long ans = f[n]; 
        long long tmp = qpow (n, n-2) * qpow (m, n-1) % mod;
        ans = (ans - tmp + mod) % mod;
        printf ("Case #%d: %lld\n", ++kase, ans);
    }
    return 0;
}