caioj1451 【快速傅里叶变换】Sigma & bzoj3527 [Zjoi2014]力

【题意】
见caioj 1451或bzoj 3527


【题解1】：推公式
快速傅立叶变换
第一步：Ei=Fi/qi，所以F计算式中的qi可以全部约掉.
第二步：简化式子，使之成为S[i]=sigma (A[j]*B[i-j])的卷积式。

设（分子：）A[i]=q[i]，（分母:）B[i]=i/(i^2)。
设（前一个累积：）C[i]=sigma(A[j]*B[i-j])，（后一个累积：）D[i]=sigma(A[n-j-1]*B[i-j])。
那么所求的E[i]=C[i]-D[i]。
不难发现C[i]已经是标准的卷积形式了，用FFT即可。
对于D[i],令A'[i]=A[n-i-1],也就是反过来，那么D[i]=sigma(A'[j]*B[i-j])，于是也用FFT即可。


【题解2】：找规律
详见http://blog.csdn.net/wzq_qwq/article/details/48155921


【小结】
学会把公式转换成卷积的形式，
当看到A[n-j-1]时，可以把A数组反过来，即A'[i]=A[n-i-1]，这是一种常规做法。

【代码】

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
 
const int maxn=1e5+10;
const double PI=acos(-1.0);
 
int R[maxn*4];
 
struct Complex
{
    double r,i;
    Complex(){}
    Complex(double _r,double _i){r=_r;i=_i;}
    friend Complex operator + (const Complex &x,const Complex &y){return Complex(x.r+y.r,x.i+y.i);}
    friend Complex operator - (const Complex &x,const Complex &y){return Complex(x.r-y.r,x.i-y.i);}
    friend Complex operator * (const Complex &x,const Complex &y){return Complex(x.r*y.r-x.i*y.i,x.r*y.i+x.i*y.r);}
}a1[maxn*4],a2[maxn*4],b[maxn*4];
 
void fft(Complex *y,int len,int on)
{
    for(int i=0;i<len;i++)
    {
        if(i<R[i]) swap(y[i],y[R[i]]);
    }
    for(int i=1;i<len;i<<=1)
    {
        Complex wn(cos(PI/i),sin(on*PI/i));
        for(int j=0;j<len;j+=i<<1)
        {
            Complex w(1,0);
            for(int k=0;k<i;k++,w=w*wn)
            {
                Complex u=y[j+k];
                Complex v=w*y[j+k+i];
                y[j+k]=u+v;
                y[j+k+i]=u-v;
            }
        }
    }
    if(on==-1)
    {
        for(int i=0;i<len;i++) y[i].r/=len;
    }
}
 
int main()
{
    int n,m,L=0;
    scanf("%d",&n);
    for(int i=0;i<n;i++)
    {
        scanf("%lf",&a1[i].r);
        a2[n-i-1]=a1[i];
    }
    for(int i=1;i<n;i++) b[i].r=1.0/i/i;
     
    m=(n-1)+(n-1);
    for(n=1;n<=m;n<<=1) L++;
    for(int i=0;i<n;i++) R[i]=(R[i>>1]>>1)|(i&1)<<(L-1);
     
    fft(a1,n,1);fft(a2,n,1);
    fft(b,n,1);
    for(int i=0;i<=n;i++)
    {
        a1[i]=a1[i]*b[i];
        a2[i]=a2[i]*b[i];
    }
    fft(a1,n,-1);fft(a2,n,-1);
     
    for(int i=0;i<=m/2;i++) printf("%.3lf\n",a1[i].r-a2[m/2-i].r);
    return 0;
}