HDU - 5730 CDQ分治 + FFT

题意：

有一串长度为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
);
}
return 0;
}