BZOJ 3527[Zjoi2014]力 FFT

题目链接:BZOJ3527
第一次学会如何写数学公式，虽然只是简单的入门，但还是有点激动。。。
首先这个题很明显是多项式乘法，但是强迫症的我过于纠结下标，以至于困惑了好久，简直SB。

注:下标均从0开始。
现在进入正题。

$$F_j=\sum_{i<j}\frac{q_i q_j} {(i-j)^2}-\sum_{i>j}\frac{q_i q_j} {(i-j)^2}$$
现在我们转化一下:记住下标从0开始
$$F_j=\sum_{i=0}^{j-1}\frac{q_i q_j}{(i-j)^2}-\sum_{i=j+1}^{n-1}\frac{q_i q_j}{(i-j)^2}$$
$$E_j=\frac{F_j}{q_j}=\sum_{i=0}^{j-1}\frac{q_i}{(i-j)^2}-\sum_{i=j+1}^{n-1}\frac{q_i}{(i-j)^2}$$

记$f(i)=q_i$，$g(i)=\frac{1}{i^2}$，其中$g(0)=0$有

$$E_j=\sum_{i=0}^{j-1}f(i)*g(j-i)-\sum_{i=j+1}^{n-1}f(i)*g(j-i)$$
等号右边第一个式子 $\sum_{i=0}^{j-1}f(i)*g(j-i)=\sum_{i=0}^jf(i)*g(j-i)$，就是卷积的形式，直接FFT算即可。

等号右边第二个式子$\sum_{i=j+1}^{n-1}f(i)*g(j-i)=\sum_{i=j}^{n-1}f(i)*g(j-i)=\sum_{i=0}^{n-j-1}f(i+j)*g(i)$
倒序处理，令$h(n-1-i-j)=f(i+j)$，有$\sum_{i=0}^{n-j-1}f(i+j)*g(i)=\sum_{i=0}^{n-j-1}h(n-1-i-j)*g(i)$，记为$X_i=\sum_{i=0}^{n-j-1}h(n-1-i-j)*g(i)$，则$X_{n-j-1}=\sum_{i=0}^{j}h(j-i)*g(i)$，也是卷积的形式，直接FFT即可。
最后$E_j=\sum_{i=0}^jf(i)*g(j-i)-X_{n-j-1}$，其中$\sum_{i=0}^jf(i)*g(j-i)$与$X_{n-j-1}$均可由FFT直接算的，题目即可解决。
代码如下:

/**************************************************************
    Problem: 3527
    User: sfailsthy
    Language: C++
    Result: Accepted
    Time:3764 ms
    Memory:17680 kb
****************************************************************/

#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdio>
#include <cmath>
using namespace std;
const int maxn=262200+10;
const double PI =acos(-1.0);

struct Complex{
    double real,image;
    Complex(){
        real=image=0.0;
    }

    Complex(double a,double b){
        real=a;
        image=b;
    }

    Complex operator + (const Complex &s) const{
        return Complex(real+s.real,image+s.image);
    }

    Complex operator - (const Complex &s) const{
        return Complex(real-s.real,image-s.image);
    }

    Complex operator * (const Complex &s) const{
        return Complex(real*s.real-image*s.image,real*s.image+image*s.real);
    }
}x1[maxn],x2[maxn],A[maxn];

int n,rev[maxn];
double a[maxn/2],b[maxn/2],c[maxn/2];

void init(int &len,int len1,int len2){
    int k=1,L=0,r,t;
    while(k<2*len1||k<2*len2){
        k<<=1;
        L++;
    }

    len=k;
    for(int i=0;i<len;i++){
        t=i,r=0,k=L;
        while(k--){
            r<<=1;
            r|=t&1;
            t>>=1;
        }
        rev[i]=r;
    }
}

void FFT(Complex x[],int len,int op){
    Complex u,t;
    for(int i=0;i<len;i++){
        A[rev[i]]=x[i];
    }
    for(int i=0;i<len;i++){
        x[i]=A[i];
    }

    for(int k=2;k<=len;k<<=1){
        Complex wn(cos(2*PI/k*op),sin(2*PI/k*op));
        for(int i=0;i<len;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;
            }
        }
    }

    if(op==-1){
        for(int i=0;i<len;i++){
            x[i].real/=len;
        }
    }
}

int main(){
    scanf("%d",&n);
    for(int i=0;i<n;i++){
        scanf("%lf",&a[i]);
    }
    for(int i=1;i<n;i++){
        b[i]=1.0/i/i;
    }

    int len,len1=n,len2=n;
    init(len,len1,len2);
    for(int i=0;i<len;i++){
        x1[i]=Complex();
        x2[i]=Complex();
    }
    for(int i=0;i<len;i++){
       x1[i]=Complex(a[i],0);
       x2[i]=Complex(b[i],0);
    }

    FFT(x1,len,1);
    FFT(x2,len,1);
    for(int i=0;i<len;i++){
        x1[i]=x1[i]*x2[i];
    }
    FFT(x1,len,-1);
    for(int i=0;i<n;i++){
        c[i]=x1[i].real;
    }

    for(int i=0;i<len;i++){
        x1[i]=Complex();
        x2[i]=Complex();
    }
    for(int i=0;i<n;i++){
        x1[i]=Complex(a[n-1-i],0);
        if(i) x2[i]=Complex(b[i],0);
    }
    FFT(x1,len,1);
    FFT(x2,len,1);
    for(int i=0;i<len;i++){
        x1[i]=x1[i]*x2[i];
    }
    FFT(x1,len,-1);
    for(int i=0;i<n;i++){
        c[i]-=x1[n-1-i].real;
    }

    for(int i=0;i<n;i++){
        printf("%.3lf\n",c[i]);
    }
    return 0;
}