bzoj 3527: [Zjoi2014]力 （FFT）

题目描述

传送门

给出$n$个数$q_i$，给出$F_j$的定义如下：
$F_j=\sum\limits_{i < j}{q_iq_j \over (i-j)^2}-\sum\limits_{i>j}{q_iq_j\over(i-j)^2}$
令$E_i=F_i/q_i$ ,求$E_i$

输入

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

输出

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

题解

FFT求卷积
$E_j=\sum\limits_{i < j}{q_i \over (i-j)^2}-\sum\limits_{i>j}{q_i\over(i-j)^2}$
我们的目标是将上面的式子化成$c_j=\sum\limits_{i=0}^j a_i*b_{j-i}$的形式，然后用傅立叶变换来求
$\sum\limits_{i < j}{q_i \over (i-j)^2}=\sum\limits_{i=0}^{j}{q_i\over(j-i)^2}$
设$f_i=q_i,g_i=1/i^2$ 其中$g_0=0$
那么上式就可以化简成$\sum\limits_{i=0}^{j}{f_i*g_{j-i}}$ ,虽然要求$i<j$,但是$g_{j-j}=0$所以这一项自动消掉了。
$\sum\limits_{i>j}{q_i\over(i-j)^2}=\sum\limits_{i=j}^n{q_i\over(i-j)^2}$
$=\sum\limits_{i=0}^{n-j}{q_{n-i}\over(j-i)^2}=\sum\limits_{i=0}^{n-j}f_{n-i}*g_{i-j}$
然后发现$n-i+i-j$是定值$n-j$,那么这一部分的答案就是FFT后$c_{n-j}$的值，但是注意这里的$f$数组是翻转后的$f$数组。

代码

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define N 300005
#define pi acos(-1)
using namespace std;
struct data{
    double x,y;
    data(double X=0,double Y=0) {
        x=X,y=Y;
    }
}f[N],f1[N],g[N],a[N],b[N];
data operator +(data a,data b){
    return data(a.x+b.x,a.y+b.y);
}
data operator -(data a,data b){
    return data(a.x-b.x,a.y-b.y);
}
data operator *(data a,data b){
    return data(a.x*b.x-a.y*b.y,a.y*b.x+a.x*b.y);
}
int n,m;
void fft(data x[N],int n,int opt)
{
    if (n==1) return;
    data l[n>>1],r[n>>1];
    for (int i=0;i<n;i+=2)
     l[i>>1]=x[i],r[i>>1]=x[i+1];
    fft(l,n>>1,opt); fft(r,n>>1,opt);
    data wn=data(cos(2*pi/n),sin(opt*2*pi/n));
    data w=data(1,0),t;
    for (int i=0;i<n>>1;i++,w=w*wn)
     t=r[i]*w,x[i]=l[i]+t,x[i+(n>>1)]=l[i]-t; 
}
int main()
{
    freopen("a.in","r",stdin);
    freopen("my.out","w",stdout);
    scanf("%d",&n); n--; m=2*n;
    for (int i=0;i<=n;i++) scanf("%lf",&f[i].x),f1[n-i].x=f[i].x;
    for (int i=1;i<=n;i++) 
     g[i].x=(double)1.0/i/i;
    for (n=1;n<=m;n<<=1);
    fft(f,n,1); fft(f1,n,1); fft(g,n,1);
    for (int i=0;i<n;i++) a[i]=f[i]*g[i];
    for (int i=0;i<n;i++) b[i]=f1[i]*g[i];
    fft(a,n,-1); fft(b,n,-1);
    for (int i=0;i<=n;i++) a[i].x=a[i].x/n;
    for (int i=0;i<=n;i++) b[i].x=b[i].x/n;
    for (int i=0;i<=m/2;i++) 
     printf("%.3lf\n",a[i].x-b[m/2-i].x);
}