[ZJOI2014] 力

第一道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;
}