HDU - 5730 CDQ分治 + FFT

最新推荐文章于 2019-08-19 15:47:30 发布

原创最新推荐文章于 2019-08-19 15:47:30 发布 · 590 阅读

4 ·

本内容遵循CC 4.0 BY-SA版权协议

标签

#acm #hdu #FFT #CDQ分治

分治同时被 2 个专栏收录

12 篇文章

订阅专栏

FFT

5 篇文章

订阅专栏

本文介绍了一种利用快速傅立叶变换(FFT)和CDQ分治技巧优化动态规划(DP)的问题解决方法。针对特定类型的DP问题，通过巧妙地结合FFT加速多项式乘法与CDQ分治策略，将复杂度降低到O(nlog^2n)，适用于求解序列划分的价值总和问题。

题意：

有一串长度为n的序列，可以任意划分，题目给出长度从1到n的块的价值，每一种划分的价值等于划分出来的所有块的价值的乘积，问所有划分方法得到的价值的总和是多少？

思路：

这题很容易就能得到状态转移方程，设dp[i]是长度为i的序列的价值总和，那么就有状态转移方程：

然而计算这样的状态转移需要花费O(n^2)的时间，不能满足要求。注意观察到这是一个离散卷积形式，可以用FFT来加速运算。

FFT可以加速多项式相乘系数的计算，这里也同样可以看作两个多项式A和B，将dp[j]看作是多项式A中j次项的系数，a[i-j]看作是多项式B中i-j次项的系数，dp[i]就可以看作是A和B相乘时得到的多项式C中i次项的系数，这样就符合了FFT的运算条件。但是如果对于每个dp[i]都用FFT运算，复杂度将达到O(n^2logn)，显然更慢，所以这里需要用到CDQ分治。

CDQ(l,r)就是用来求出l到r的dp值，每次二分区间，先递归计算出(l,m)的dp值，再利用(l,m)的dp值来更新(m+1,r)的dp值，更新的过程使用FFT加速运算，算法总复杂度为O(nlognlogn)。

这里在计算(l,m)的贡献的时候，因为要使用FFT，需要将下标统一到1到len存在x1和x2中，注意细节。

代码：

#include <bits/stdc++.h>
using namespace std;
const double pi = acos(-1);
const int MAXN = 555555;
const int MOD = 313;

struct plex {       // 定义复数类
    double x, y;
    plex (double _x = 0.0, 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);
    }
};

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

void fft(plex y[], int len, int on) {   // FFT过程，on==1时，将系数表达转换成点值表达，on==-1时，将点值表达转换成系数表达
    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;
    }
}

int a[MAXN], dp[MAXN];
plex x1[MAXN], x2[MAXN];

void cdq(int l, int r) {
    if (l == r) return;
    int m = (r - l) / 2 + l;
    cdq(l, m);
    int len = 2;
    while (len < (r - l + 1) * 2) len <<= 1;
    for (int i = 0; i < len; i++)
        x1[i] = x2[i] = plex(0, 0);
    for (int i = l; i <= m; i++)
        x1[i - l + 1] = plex(dp[i], 0);
    for (int i = 1; i <= r - l; i++)
        x2[i] = plex(a[i], 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 = m + 1; i <= r; i++)
        dp[i] = (dp[i] + (int)(x1[i - l + 1].x + 0.5)) % MOD;
    cdq(m + 1, r);
}

int main () {
    int n;
    while (scanf("%d", &n), n) {
        memset(dp, 0, sizeof(dp));
        for (int i = 1; i <= n; i++) {
            scanf("%d", &a[i]);
            a[i] %= MOD;
            dp[i] = a[i];
        }
        cdq(1, n);
        printf("%d\n", dp[n]);
    }
    return 0;
}