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学习资料:
从多项式乘法到快速傅里叶变换
<从入门到崩溃>快速傅里叶变换学习笔记


代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
/*************************************************************************
> 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;
}
文章目录
  1. 1. 题目链接:
  2. 2. 题意:
  3. 3. 分析:
  4. 4. 代码: