BZOJ3527ZJOI力——FFT

Description

给出n个数qi，给出Fj的定义如下：

令Ei=Fi/qi，求Ei.

Input
第一行一个整数n。
接下来n行每行输入一个数，第i行表示qi。
n≤100000，0< qi< 1000000000

Output
n行，第i行输出Ei。与标准答案误差不超过1e-2即可。

Sample Input
5
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880

Sample Output
-16838672.693
3439.793
7509018.566
4595686.886
10903040.872

---

这道题可谓是FFT的模板题，我们可以化一下式子：

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

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

我们令$g[i]=\frac{1}{i^2}$,则$\sum_{i<j}\frac{q_i}{(i-j)^2}=\sum_{i<j}q_i*g_{j-i}=\sum_{i=0}^{j-1}q_i*g_{j-i}$

我们发现$\sum_{i=0}^{j-1}q_i*g_{j-i}$就是卷积的形式，所以我们可以用FFT来处理。

现在的关键就是处理$\sum_{i>j}\frac{q_i}{(i-j)^2}$

$\sum_{i>j}\frac{q_i}{(i-j)^2}=\sum_{i=j+1}^{n-1}q_i*g(j-i)=\sum_{i=j}^{n-1}q_i*g(j-i)=\sum_{i=0}^{n-j-1}q_{i+j}*g_{i}$我们要想办法让其变为卷积的形式，我们设$h_{n-i-j-1}=q_{i+j}$，则原式等于$\sum_{i=0}^{n-j-1}h_{n-1-i-j}*g{i}$，计作$X_i$,那么$X_{n-j-1}=h_{j-i}*g_{i}$也满足卷积的形式，可以用FFT处理。

所以我们用FFT处理出两个卷积，然后求出答案即可。

#include<bits/stdc++.h>
#define db double
#define MAXN 300000
using namespace std;
int read(){
    char c;int x=0,y=1;while(c=getchar(),(c<'0'||c>'9')&&c!='-');
    if(c=='-') y=-1;else x=c-'0';while(c=getchar(),c>='0'&&c<='9')
    x=x*10+c-'0';return x*y;
}
const db pi=acos(-1.0);
struct Complex{
    db x,y;
    void clear(){x=0,y=0;}
    Complex (db xx=0,db yy=0){x=xx,y=yy;}
}F[MAXN],H[MAXN];
Complex operator+(Complex a,Complex b){return Complex(a.x+b.x,a.y+b.y);}
Complex operator-(Complex a,Complex b){return Complex(a.x-b.x,a.y-b.y);}
Complex operator*(Complex a,Complex b){return Complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
int n,limit=1,l,r[MAXN];
db a[MAXN],b[MAXN],c[MAXN];
void FFT(Complex *A,int type){
    for(int i=0;i<limit;i++) if(i<r[i]) swap(A[i],A[r[i]]);
    for(int mid=1;mid<limit;mid<<=1){
        Complex Wn(cos(pi/mid),type*sin(pi/mid));
        for(int R=mid<<1,j=0;j<limit;j+=R){
            Complex w(1,0);
            for(int k=0;k<mid;k++,w=w*Wn){
                Complex x=A[j+k],y=w*A[j+k+mid];
                A[j+k]=x+y;A[j+k+mid]=x-y;
            }
        }
    }
}
int main()
{
    n=read();
    for(int i=0;i<n;i++) scanf("%lf",&a[i]);
    for(int i=1;i<n;i++) b[i]=(db)1/i/i;  //写成1/(i*i)会爆int
    while(limit<=n+n-2) limit<<=1,l++;
    for(int i=0;i<limit;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    for(int i=0;i<limit;i++) F[i].clear(),H[i].clear();
    for(int i=0;i<limit;i++) F[i].x=a[i],H[i].x=b[i];
    FFT(F,1);FFT(H,1);
    for(int i=0;i<limit;i++) F[i]=F[i]*H[i];
    FFT(F,-1);
    for(int i=0;i<n;i++) c[i]=F[i].x/limit+0.5;
    for(int i=0;i<limit;i++) F[i].clear(),H[i].clear();
    for(int i=0;i<n;i++){
        H[i].x=b[i],0;F[i].x=a[n-1-i],0;
    }
    FFT(F,1);FFT(H,1);
    for(int i=0;i<limit;i++) F[i]=F[i]*H[i];
    FFT(F,-1);
    for(int i=0;i<n;i++){
        c[i]-=F[n-1-i].x/limit+0.5;
    }
    for(int i=0;i<n;i++) printf("%.3lf\n",c[i]);
    return 0;
}