第一道FFT例题,代码和脑子都一片混乱,感觉难想,细节很多
传送门:$>here<$
题意
求$E_k=\sum\limits_{i<k}\dfrac{q_i}{(i-k)^2}-\sum\limits_{i>k}\dfrac{q_i}{(i-k)^2}$
数据范围:$n \leq 1e5$
Solution
之前我们用FFT完成多项式乘法时,其实就是求两个系数向量的卷积。$A$与$B$乘积的系数向量具有如下的卷积形式:$$C_k=\sum\limits_{i \leq k}a_ib_{k-i}$$于是我们希望把题目里的这个式子拆成卷积的形式,这样就可以FFT了。
我们先用系数来表示,令$f(x)=\frac{1}{x^2}$。设$P_k=\sum\limits_{i<k}q(i)f(k-i)$,$Q_k=\sum\limits_{i>k}q(i)f(k-i)$,则$E_k=P_k-Q_k$
由于$f(0)$不存在,我们可以设$f(0)=0$,这样不会影响结果。这样的话$P$和$Q$的范围都可以取等号了。
$P_k=\sum\limits_{i \leq k}q(i)f(k-i)$,$Q_k=\sum\limits_{i \geq k}q(i)f(k-i)$
$P_k$就是标准的卷积形式,因此直接由$q(i)=q_i$,$f(i)=\frac{1}{i^2}$作为系数向量FFT一遍求出。
难点在于$Q_k$的求解。卷积形式的范围是从小到大的。因此可以换元,令$i=n-j$,代入得
$Q_k=\sum\limits_{j \leq n-k}q(n-j)f(k-n+j)$
这样一来不能肯定不能再用$q$或者$f$这两个函数了,因为它们的下标之和应当等于$n-k$
由于$f$是偶函数,所以$f(k-n+j)=f(n-j-k)$。如果$q$的下标是$j$就行了。于是设$q'(j)=q(n-j)$,就有
$Q_k=\sum\limits_{j \leq n-k}q'(j)f(n-j-k)$
因此我们所需要做的就是将原$q$数组镜像(倒序),然后再做一次卷积。$Q_k$就是下标为$n-k$的系数。
my code
注意下标从0开始。
/*By DennyQi 2018*/
#include <cstdio>
#include <queue>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
typedef long long ll;
const int MAXN = 600010;
const int INF = 0x3f3f3f3f;
const double Pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164;
inline int Max(const int a, const int b){ return (a > b) ? a : b; }
inline int Min(const int a, const int b){ return (a < b) ? a : b; }
inline int read(){
int x = 0; int w = 1; register char c = getchar();
for(; c ^ '-' && (c < '0' || c > '9'); c = getchar());
if(c == '-') w = -1, c = getchar();
for(; c >= '0' && c <= '9'; c = getchar()) x = (x<<3) + (x<<1) + c - '0'; return x * w;
}
struct Complex{
double a,b;
Complex(){
a = 0.0, b = 0.0;
}
Complex(const double x, const double y){
a = x, b = y;
}
Complex operator + (const Complex& u){
return Complex(a+u.a,b+u.b);
}
Complex operator - (const Complex& u){
return Complex(a-u.a,b-u.b);
}
Complex operator * (const Complex& u){
return Complex(a*u.a-b*u.b,a*u.b+b*u.a);
}
};
int N,Lim,lg;
int r[MAXN];
double q[MAXN],ans1[MAXN],ans2[MAXN];
Complex a[MAXN],b[MAXN];
inline void FFT(int n, Complex* a, int cvs){
for(int i = 0; i < n; ++i){
if(i < r[i]) swap(a[i],a[r[i]]);
}
Complex A1,A2,Wn,w;
for(int i = 1; i <= n; i<<=1){
for(int j = 0; j < n; j+=i){
Wn=Complex(1,0),w=Complex(cos(2*Pi/i*cvs),sin(2*Pi/i*cvs));
for(int k = 0; k < (i/2); ++k){
A1 = a[j+k], A2 = a[j+k+(i/2)];
a[j+k] = A1 + Wn*A2;
a[j+k+(i/2)] = A1 - Wn*A2;
Wn = Wn * w;
}
}
}
}
int main(){
// freopen(".in","r",stdin);
scanf("%d",&N);
for(int i = 1; i <= N; ++i){
scanf("%lf",q+i);
}
Lim = 1;
while(Lim <= (N<<1)) Lim<<=1, ++lg;
for(int i = 0; i < Lim; ++i){
r[i] = (r[i>>1]>>1) | ((i&1)<<(lg-1));
}
for(int i = 0; i < Lim; ++i){
if(i <= N && i > 0){
a[i].a = q[i];
b[i].a = (double)(1 / ((double)(i)*(double)(i)));
}
}
FFT(Lim,a,1);
FFT(Lim,b,1);
for(int i = 0; i < Lim; ++i){
a[i] = a[i] * b[i];
}
FFT(Lim,a,-1);
for(int i = 0; i < Lim; ++i){
ans1[i] = a[i].a / Lim;
}
for(int i = 0; i < Lim; ++i){
a[i] = Complex(0,0);
b[i] = Complex(0,0);
if(N-i>=0){
a[i].a = q[N-i];
}
if(i <= N && i > 0){
b[i].a = (double)(1 / ((double)(i)*(double)(i)));
}
}
FFT(Lim,a,1);
FFT(Lim,b,1);
for(int i = 0; i < Lim; ++i){
a[i] = a[i] * b[i];
}
FFT(Lim,a,-1);
for(int i = 0; i < Lim; ++i){
ans2[i] = a[i].a / Lim;
}
for(int i = 1; i <= N; ++i){
printf("%.3f\n",ans1[i]-ans2[N-i]);
}
return 0;
}