[ZJOI2014]力 FFT基本套路实践

此文中介绍了FFT的基本套路：

[BZOJ]4503 两个串：我的第一次FFT尝试

本题便是对FFT基本套路的又一次实践，即翻转一个串，把对应相乘被转化成卷积。

题目描述

给出$n$个数$q_i$，$F_j$的定义如下：

$$F_j=\sum_{i<j}{\frac{q_iq_j}{(i-j)^2}}-\sum_{i>j}{\frac{q_iq_j}{(i-j)^2}}$$

令：$E_i=F_i/q_i$，求$E_i$。

输入格式

第一行一个数n，接下来n行每行一个数，第i行的数为$q_i$。

输出格式

n行，第i行的数表示$E_i$。

样例输入

5
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880

样例输出

-16838672.693
3439.793
7509018.566
4595686.886
10903040.872

解题思路

把$F_j$中的所有$q_j$都提出来，有：

$$F_j=q_j \cdot (\sum_{i<j}{\frac{q_i}{(i-j)^2}}-\sum_{i>j}{\frac{q_i}{(i-j)^2}})$$

所以有:

$$E_i=\frac{F_i}{q_i}=\sum_{i<j}{\frac{q_i}{(i-j)^2}}-\sum_{i>j}{\frac{q_i}{(i-j)^2}}$$

把这个式子拆成两半，令：

$$f(j)=\sum_{i<j}{\frac{q_i}{(i-j)^2}}=\sum_{i=0}^{j-1}{\frac{q_i}{(i-j)^2}}$$
$$g(j)=\sum_{i>j}{\frac{q_i}{(i-j)^2}}=\sum_{i=j+1}^{n-1}{\frac{q_i}{(i-j)^2}}=\sum_{i=j+1}^{n}{\frac{q_i}{(i-j)^2}}$$

则有:

$$E_i=f(i)-g(i)$$

先看$f(i)$：$f(j)=\sum_{i=0}^{j-1}{\frac{q_i}{(i-j)^2}}$，不妨令$i=j-i$（调整循环顺序，结果不变），有：

$$f(j)=\sum_{(j-i)=0}^{(j-i)=j-1}{\frac{q_{j-i}}{((j-i)-j)^2}}=\sum_{i=1}^{j}{\frac{q_{j-i}}{i^2}}$$

定义向量A，B，使得：$A_i=q_i$，$B_i=\frac{1}{i^2}(i \in 1..n)$，特别地，$B_0=0$。则有：

$$f(j)=\sum_{i=1}^{j}{q_{j-i}\cdot \frac{1}{i^2}}=\sum_{i=0}^{j}{A_{j-i}\cdot B_{i}}$$

令向量C为，向量A、B的卷积。即：$C=A \times B$，有：

$$C_k=\sum_{i+j=k}{A_j \cdot B_i}=\sum_{i=0}^{k}{A_{k-i} \cdot B_i}$$

由此证明：

$$f(j)=C_j=(A \times B)_j$$

然后考虑$g(j)$：$g(j)=\sum_{i=j+1}^{n}{\frac{q_i}{(i-j)^2}}$不妨令$i={j-i}$（同上文，调整循环顺序，结果不变），有：

$$g(j)=\sum_{{j-i}=j+1}^{{j-i}=n}{\frac{q_{j-i}}{(({j-i})-j)^2}}=\sum_{i=j-n}^{-1}{\frac{q_{j-i}}{(-i)^2}}$$

i小于零看着太别扭，再令$i=-i$，有：

$$g(j)=\sum_{i=1}^{n-j}{\frac{q_{j+i}}{i^2}}=\sum_{i=1}^{n-j}{q_{j+i} \cdot \frac{1}{i^2}}=\sum_{i=0}^{n-j}{A_{j+i} \cdot B_i}$$

此时，我们发现，i与j+i的差是一个定值。

定义一个向量$B'$，使得$B'_i=B_{n-i}$，即对B向量进行翻转。那么就有：

$$g(j)=\sum_{i=0}^{n-j}{A_{j+i} \cdot B'_{n-i}}$$

这次，不难看出，n-i与j+i的和是一个定值了。

定义向量$C'$，有$C'=A \times B'$，即：

$$C'_k=\sum _{i+j=k}{A_i \cdot B'_j}=\sum_{i=0}^{k}{A_{i}\cdot B'_{k-i}}$$

那么：

$$C'_{n+j}=\sum_{i=0}^{n+j}{A_{i}\cdot B'_{n+j-i}}$$

令$i=j+i$，有：

$$C'_{n+j}=\sum_{j+i=0}^{j+i=n+j}{A_{j+i}\cdot B'_{n+j-({j+i})}}=\sum_{i=-j}^{n}{A_{j+i}\cdot B'_{n-i}}$$

当$i<0$时，$n-i>n$，$B'_{n-i}=0$，可以不考虑；当$i>n-j$时，$j+i>n$，$A_{j+i}=0$，也可以不考虑，所以：

$$C'_{n+j}=\sum_{i=-j}^{n}{A_{j+i}\cdot B'_{n-i}}=\sum_{i=0}^{n-j}{A_{j+i}\cdot B'_{n-i}}=g(j)$$

结合上文中所述，我们可以利用FFT以$O(nlgn)$的时间复杂度计算卷积$C=A \times B$和$C'=A \times B'$ 。输出的结果就是$E_i=C_{i}-C'_{n+i}$ 。

代码

#include<algorithm>
#include<cmath>
#include<complex>
using namespace std;
typedef complex<double>cd;

#include<stdio.h>
#include<stdlib.h>

#define maxl ((100000+3)*4)
cd A[maxl],B[maxl],invB[maxl];
cd f[maxl],g[maxl];

double ans[maxl];

int rev[maxl];
void get_rev(int bit){
    for(int i=0;i<(1<<bit);i++){
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    }
}

void fft(cd* a,int n,int dft){
    for(int i=0;i<n;i++){
        if(i<rev[i])
            swap(a[i],a[rev[i]]);
    }
    for(int step=1;step<n;step<<=1){
        cd wn=exp(cd(0,dft*M_PI/step));
        for(int j=0;j<n;j+=step<<1){
            cd wnk(1,0);
            for(int k=j;k<j+step;k++){
                cd x=a[k];
                cd y=wnk*a[k+step];
                a[k]=x+y;
                a[k+step]=x-y;
                wnk*=wn;
            }
        }
    }
    if(dft==-1){
        for(int i=0;i<n;i++)
            a[i]/=n;
    }
}

int main(){
    int n;
    scanf("%d",&n);
    for(int i=0;i<n;i++){
        double x;
        scanf("%lf",&x);
        A[i]=x;
    }
    for(int i=1;i<=n;i++)B[i]=1.0/i/i;
    for(int i=0;i<=n;i++)invB[i]=B[n-i];
    int bit,s;
    for(bit=1,s=2;(1<<bit)<(n+(n+1))-1;bit++)s<<=1;
    get_rev(bit);
    fft(A,s,1);fft(B,s,1);fft(invB,s,1);
    for(int i=0;i<s;i++){
        f[i]=A[i]*B[i];
        //printf("f[%d]=(%lf,%lf)\n",i,f[i].real(),f[i].imag());
        g[i]=A[i]*invB[i];
    }
    fft(f,s,-1);fft(g,s,-1);
    for(int i=0;i<n;i++){
        ans[i]=f[i].real()-g[n+i].real();
        printf("%.3lf\n",ans[i]);
    }
    return 0;
}