关键词:快速傅里叶变换(FFT)  卷积 相关题目:[ZJOI2014] 力

「ZJOI2014」力

Time Limit: 3Sec   Memory Limit: 256MB   Special judge

题目描述

给出 nn 个数 qiq_i,给出 FjF_j 的定义如下:

Fj=i<jqiqj(ij)2i>jqiqj(ij)2F_{j}=\sum_{i<j}\frac{q_{i}q_{j}}{(i-j)^{2}}-\sum_{i>j}\frac{q_{i}q_{j}}{(i-j)^{2}}

Ei=Fi/qiE_{i}=F_{i}/q_{i} ,求 EiE_{i}

输入输出格式

输入格式: 第一行一个整数nn。 接下来 nn 行每行输入一个数,第 ii 行表示 qiq_i

输出格式: nn 行,第 ii 行输出 EiE_i。与标准答案误差不超过 10210^{-2} 即可。

分析及解决

化简 EjE_{j}

Ej=Fj/qj=i<jqi(ij)2i>jqi(ij)2E_{j}=F_{j}/q_{j}=\sum_{i<j}\frac{q_{i}}{(i-j)^{2}}-\sum_{i>j}\frac{q_{i}}{(i-j)^{2}}

fi=qif_{i}=q_{i}gi=1i2g_{i}=\frac{1}{i^{2}} ,特别地,约定 g0=0g_{0}=0

对于左边的式子:

这一部分我们可以很容易的通过 快速傅里叶变换 求得。

对于右边的式子:

注意到 njn-j 为定值,于是这一部分同样可以通过 快速傅里叶变换 求得,这里要注意这里的 ff 为翻转后的 ff

两遍FFT,完了......

1
2
**题卡精度 (差评
特别注意这一段: cpx(1.0/(double)i/(double)i, 0.0);

AC代码

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
#include <cstdio>
#include <cstring>
#include <cmath>
#include <cctype>
#include <memory>
#define re register int
const int N=2097152+5, M=1000020;
const double pi=acos(-1);
struct cpx {
double r,i;
cpx(){}
cpx(double _r,double _i){r=_r, i=_i;}
inline cpx operator * (const cpx& rhy)const{
return cpx(r*rhy.r - i*rhy.i, r*rhy.i + i*rhy.r);
}
inline cpx operator *= (const cpx& rhy){
*this = *this * rhy;
}
inline cpx operator + (const cpx& rhy)const{
return cpx(r+rhy.r, i+rhy.i);
}
inline cpx operator - (const cpx& rhy)const{
return cpx(r-rhy.r, i-rhy.i);
}
}f[N], g[N], c[N]; //complex

int n, R[N];
double q[M], ans[M];

inline void FFT(int n, cpx *a, int f){
for(re i = 0 ; i < n ; ++i)
if(R[i] > i) std::swap(a[R[i]], a[i]); //二进制反转

for(re l = 1 ; l < n ; l<<=1) {
cpx wn(cos(pi/l), f*sin(pi/l));
for(re j = 0 ; j < n ; j += (l<<1)) {
cpx omega(1.0, 0.0);
for(re k = 0 ; k < l ; ++k,omega*=wn) {
cpx t = omega * a[j+k+l];
a[j+k+l] = a[j+k] - t;
a[j+k] = a[j+k] + t;
}
}
} // 分治+蝴蝶变换

if(f == -1)
for(re i = 0 ; i < n ; ++i) a[i].r /= (double)n; //系数转换
}

int main(int argc,char *argv[]){
scanf("%d", &n);
int x;
for(re i = 0 ; i < n ; ++i) scanf("%lf", &q[i]);

int _n=1;
while(_n < n*2) _n<<=1;
for(re i = 0 ; i < _n ; ++i)
R[i]=( (R[i>>1]>>1) | (_n>>(i&1)) ) & (_n-1); //二进制反转预处理

for(re i = 0 ; i < n ; ++i) f[i] = cpx(q[i], 0.0);
for(re i = 1 ; i < n ; ++i) g[i] = cpx(1.0/(double)i/(double)i, 0.0);
FFT(_n, f, 1), FFT(_n, g, 1);
for(re i = 0 ; i < _n ; ++i) c[i] = f[i]*g[i]; //点值转换
FFT(_n, c, -1);
for(re i = 0 ; i < n ; ++i) ans[i] = c[i].r;

memset(f, 0, sizeof f), memset(g, 0, sizeof g);

for(re i = 0 ; i < n ; ++i) f[i] = cpx(q[n-i-1], 0.0);
for(re i = 1 ; i < n ; ++i) g[i] = cpx(1.0/(double)i/(double)i, 0.0);
FFT(_n, f, 1), FFT(_n, g, 1);
for(re i = 0 ; i < _n ; ++i) c[i] = f[i]*g[i];
FFT(_n, c, -1);
for(re i = 0 ; i < n ; ++i) ans[i] -= c[n-i-1].r;

for(re i = 0 ; i < n ; ++i) printf("%.3f\n", ans[i]);
}

或者另一种写法

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
#include <cstdio>
#include <cmath>
#include <memory>
#define re register int
const int N=2097152+5, M=1000020;
const double pi=acos(-1);
struct cpx {
double r,i;
cpx(){}
cpx(double _r,double _i){r=_r, i=_i;}
inline cpx operator * (const cpx& rhy)const{
return cpx(r*rhy.r - i*rhy.i, r*rhy.i + i*rhy.r);
}
inline cpx operator *= (const cpx& rhy){
*this = *this * rhy;
}
inline cpx operator + (const cpx& rhy)const{
return cpx(r+rhy.r, i+rhy.i);
}
inline cpx operator - (const cpx& rhy)const{
return cpx(r-rhy.r, i-rhy.i);
}
}a[N], b[N];

int n, R[N];
double c[M];

inline void FFT(int n, cpx *a, int f){
for(re i = 0 ; i < n ; ++i)
if(R[i] > i) std::swap(a[R[i]], a[i]);

for(re l = 1 ; l < n ; l<<=1) {
cpx wn(cos(pi/l), f*sin(pi/l));
for(re j = 0 ; j < n ; j += (l<<1)) {
cpx omega(1.0, 0.0);
for(re k = 0 ; k < l ; ++k,omega*=wn) {
cpx t = omega * a[j+k+l];
a[j+k+l] = a[j+k] - t;
a[j+k] = a[j+k] + t;
}
}
}

if(f == -1)
for(re i = 0 ; i < n ; ++i) a[i].r /= (double)n;
}

int main(int argc,char *argv[]){
scanf("%d", &n);
int x;
for(re i = 0 ; i < n ; ++i) scanf("%lf", &a[i].r);
for(re i = 0 ; i < n-1 ; ++i)
b[i].r = -1.0/(double)(n-i-1)/(double)(n-i-1), b[(n-1)*2-i].r = -b[i].r;
int n_ = n; n = (n-1)*2+1;

int _n=1;
while(_n < n*2) _n<<=1;
for(re i = 0 ; i < _n ; ++i)
R[i]=( (R[i>>1]>>1) | (_n>>(i&1)) ) & (_n-1);

FFT(_n, a, 1), FFT(_n, b, 1);
for(re i = 0 ; i < _n ; ++i) a[i]*=b[i];
FFT(_n, a, -1);
for(re i = 0 ; i < n ; ++i) c[i] = a[i].r;

for(re i = n_-1 ; i < n ; ++i) printf("%.3f\n", c[i]);
}