2016-08-24

BZOJ 3527 [ZJOI2014]力【FFT】

题目链接：

http://www.lydsy.com/JudgeOnline/problem.php?id=3527

题意：

BZOJ上没有题面，所以在这粘一下。

Description
给出$n$个数$q_i$，给出$Fj$的定义如下：
$F_j=\sum_{i\lt j}\frac{q_iq_j}{ (i-j)^{2} }-\sum_{i\gt j}\frac{q_iq_j}{ (i-j)^{2} }$
令$Ei=\frac{Fi}{q_i}$，求$Ei$.
Input
第一行一个整数$n$。
接下来$n$行每行输入一个数，第$i$行表示$q_i$。
Output
$n$行，第$i$行输出$E_i$。
与标准答案误差不超过$1e-2$即可。
Sample Input
5
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880
Sample Output
-16838672.693
3439.793
7509018.566
4595686.886
10903040.872
Hint
$对于30\%的数据，n≤1000$。
$对于50\%的数据，n≤60000$。
$对于100\%的数据，n≤100000，0<q_i<1000000000$。
Source
感谢nodgd放题

分析：

化简成比较直观的形式就是
$$ E_j=\sum_{i=1}^{j-1}\frac{q_{i}}{(i-j)^{2}}-\sum_{i=j+1}^{n}\frac{q_i}{(i-j)^{2}}$$
令$f(i) = q_i,\ g(i) = \frac{1}{i^2}$，
则$$E_j=\sum_{i=1}^{j-1}f(i) \times g(j-i)-\sum_{i=j+1}^nf(i) \times g(j-i)$$
这样我们发现左右两边都是很卷积的形式，$c_i = \sum_{j=0}^ia_j\times b_{i-j}$，前一部分直接$FFT$算。
后一部分，$\sum_{i=j+1}^nf(i) \times g(j-i)=\sum_{i=1}^{n-j}f(i+j) \times g(i)$
倒序处理并使下标从$0$开始，令$f’(n-i-j)=f(i+j)$，则第二部分变为$\sum_{i=0}^{n-j-1}f’(n-i-j-1) \times g(i)$，转化为卷积的形式用$FFT$解即可。

非常好的FFT学习资料：
从多项式乘法到快速傅里叶变换
 <从入门到崩溃>快速傅里叶变换学习笔记

代码：

/*************************************************************************
    > File Name: 3527.cpp
    > Author: jiangyuzhu
    > Mail: 834138558@qq.com 
    > Created Time: 2016/8/24 11:24:18
 ************************************************************************/
#include<iostream>
#include<cstdio>
#include<cstring>
#include<queue>
#include<vector>
#include<set>
#include<map>
#include<cmath>
#include<algorithm>
using namespace std;
const int maxn = 4e5 + 5;
const double pi = acos(-1.0); 
struct complex
{
    double r, i;
    complex() {r = i = 0;}
    complex(double a,double b){r = a; i = b;}
    complex operator + (const complex &o) const{
        return complex(r + o.r, i + o.i);
    }
    complex operator - (const complex &o) const{
        return complex(r - o.r, i - o.i);
    }
    complex operator * (const complex &o) const{
        return complex(r * o.r - i * o.i, r * o.i + i * o.r);
    }
}x1[maxn], x2[maxn], A[maxn];
int n, rev[maxn];
double a[maxn], b[maxn], c[maxn];
void init(int &l)
{
	int L, k, r, t;
	k = 1; L = 0;
    while(k < l) k <<= 1, ++ L;
    l = k;
    for(int i = 0; i < l; i++){
        t = i; r = 0; k = L;
        while(k--){
            r <<= 1;
            r |= t & 1; 
            t >>= 1;
        }
        rev[i] = r;
    }
}
void fft(complex *x, int l, int op)
{
    complex u, t;
    for(int i = 0; i < l; i++) A[rev[i]] = x[i];
    for(int i = 0; i < l; i++) x[i] = A[i];
    for(int k = 2; k <= l; k <<= 1){
        complex wn(cos(2 * pi/ k * op), sin(2 * pi / k * op));
        for(int i = 0; i < l; i += k){
            complex w(1, 0);
            for(int j = 0; j < k/2; j++){
                u = x[i + j]; t = w * x[i + j + k/2];
                x[i + j] = u + t; x[i + j + k/2] = u - t;
                w = w * wn;
            }
        }
    }
	//IDFT
    for(int i = 0; i < l && op == -1; i++)  x[i].r /= l;
}
int main(void)
{
	int l;
    scanf("%d", &n); l = n * 2 - 1;
    for(int i = 0; i < n; i++) 	scanf("%lf", &a[i]);
    init(l);
	for(int i = 0; i < n; i++) x1[i] = complex(a[i], 0);
    for(int i = 1; i < n; i++){
		b[i] = 1.0 / i / i;
		x2[i] = complex(b[i], 0);
	}
    fft(x1, l, 1); fft(x2, l, 1);
    for(int i = 0; i < l; i++)	x1[i] = x1[i] * x2[i];	
    fft(x1, l, -1);
    for(int i = 0; i < n; i++){
		c[i] = x1[i].r;
	}
	for(int i = 0; i < l; i++){
		x1[i] = complex(0,0);
		x2[i] = complex(0,0);
	}
	for(int i = 0; i < n; i++){
		x1[i] = complex(a[n - i - 1], 0);
	}
	for(int i = 1; i < n; i++){
		x2[i] = complex(b[i], 0);
	}
	fft(x1, l, 1); fft(x2, l, 1);
    for(int i = 0; i < l; i++)	x1[i] = x1[i] * x2[i];
    fft(x1, l, -1);
    for(int i = 0; i < n; i++) c[i] -= x1[n - i - 1].r;
	for(int i = 0; i < n; i++){
		printf("%.3f\n", c[i]);
	}
    return 0;
}