[BZOJ3527][Zjoi2014]力（FFT）

=== ===

这里放传送门

=== ===

题解

一开始脑子叉劈了，全程思考怎么把这两个求和的东西合并起来变成一个卷积，然后就GG了= =但是如果先把$E_i$的柿子给画出来：

$$E_i=\sum_{j \lt i}{q_j \over {(i-j)^2}}-\sum_{j \gt i}{q_j \over {(i-j)^2}}$$
可以发现这两大块求和的式子都长得好像一个卷积。。然而第二块因为是从大到小循环而不是从小到大循环看起来好像比较GG。。那么我们可以把它倒过来，也就是构造一个新的 $q'_i$，使得 $q'_i=q_{n-i}$，然后我们这个式子就可以变成： $$E_i=\sum_{j=0}^{i-1}{q_j \over {(i-j)^2}}-\sum_{j=0}^{n-i-1}{q'_j \over {(n-i-j)^2}}$$
设 $f(i)=q_i$， $g(i)={1 \over i^2}$， $f'(i)=q'_i$，上面那个式子就可以非常愉快地写成： $$E_i=\sum_{j=0}^{i-1}f(j)g(i-j)-\sum_{j=0}^{n-i-1}f'(j)g(n-i-j)$$
可是问题又来了，卷积的话每一个求和函数都是要从0枚举到i的，也就是卷积的形式应该是 $h(i)=\sum\limits_{j=0}^if(j)g(i-j)$，但是这里因为 $g$函数在自变量为0的位置没有定义，所以我们枚举的时候只枚举到了i-1。不过解决的办法就是直接令$g(0)=0$，那么多枚举的那一项对最后结果肯定是没有贡献的。两次FFT就搞出来了。

那个。。还有，ATP一开始写的时候WA了一坨。。搞到数据以后发现最终结果比正确答案差了好几十。。中间结果能差到好几十万。。N久查错未果后懵逼的ATP把计算g函数的时候写的

for (int i=1;i<=n;i++) b[i].x=1/(double)(i*i);

改成了

for (int i=1;i<=n;i++) b[i].x=1/((double)i*(double)i);

然后就对了？！！就TMD对了？！！！！What the Fuck？！！！

代码

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
const double PI=4*atan(1);
int n,R[300010];
double h1[300010],h2[300010],q[300010];
struct Complex{
    double x,y;
    Complex(double X=0,double Y=0){x=X;y=Y;}
    Complex operator + (const Complex &a){return Complex(x+a.x,y+a.y);}
    Complex operator - (const Complex &a){return Complex(x-a.x,y-a.y);}
    Complex operator * (const Complex &a){return Complex(x*a.x-y*a.y,x*a.y+y*a.x);}
}a[300010],b[300010];
void FFT(Complex *a,int N,int opt){
    Complex w,wn,x,y;
    for (int i=0;i<N;i++)
      if (i<R[i]) swap(a[i],a[R[i]]);
    for (int k=1;k<N;k<<=1){
        wn=Complex(cos(PI/k),opt*sin(PI/k));
        for (int p=(k<<1),i=0;i<N;i+=p){
            w=Complex(1,0);
            for (int j=0;j<k;j++){
                x=a[i+j];y=w*a[i+j+k];
                a[i+j]=x+y;a[i+j+k]=x-y;
                w=w*wn;
            }
        }
    }
}
void Calc(double *h,int n){ 
    int m=n+n,L=0;
    for (n=1;n<=m;n<<=1) L++;
    memset(R,0,sizeof(R));
    for (int i=0;i<n;i++)
      R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
    FFT(a,n,1);FFT(b,n,1);
    for (int i=0;i<=n;i++) a[i]=a[i]*b[i];
    FFT(a,n,-1);
    for (int i=0;i<=n;i++) h[i]=a[i].x/(double)n;
}
int main()
{
    scanf("%d",&n);
    for (int i=1;i<=n;i++) scanf("%lf",&q[i]);
    a[0]=b[0]=Complex();
    for (int i=1;i<=n;i++) a[i].x=q[i];
    for (int i=1;i<=n;i++) b[i].x=1/((double)i*(double)i);
    Calc(h1,n);
    memset(a,0,sizeof(a));memset(b,0,sizeof(b));
    for (int i=0;i<=n;i++) a[i]=Complex(q[n-i],0);
    for (int i=1;i<=n;i++) b[i].x=1/((double)i*(double)i);
    Calc(h2,n);
    for (int i=1;i<=n;i++)
      printf("%lf\n",h1[i]-h2[n-i]);
    return 0;
}

偏偏在最后出现的补充说明

MD以后强转老老实实在每一个变量前都加括号= =
FFT这玩意儿确定不是在考脑洞吗。。。