BZOJ3527 力

题目

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

令Ei=Fi/qi，求Ei.

思路

年轻人的第一道FFT。
\[ E_j=\sum_{i<j}\frac{q_i}{(i-j)^2}-\sum_{i>j}\frac{q_i}{(i-j)^2} \]
设：
\[ f[i]=q_i\\ g[i]=\frac{1}{i^2} \]
那么原式可以化为：
\[ E_j=\sum_{i<j}f[i]g[j-i]-\sum_{i>j}f[i]g[i-j] \]
也就是:
\[ E_j=\sum_{i=0}^{j-1}f[i]g[j-i]-\sum_{i=j+1}^{n}f[i]g[i-j] \]
前半部分是卷积的形式，显然可以用FFT维护，重点是后面的部分。

我们来变一下形：
\[ \sum_{i=j+1}^{n}f[i]g[i-j]\\ 令i=i-j\\ \sum_{i=0}^{n-j}f[i+j]g[i] \]
然后将\(i+j+i\)凑成\(n-j\),也就是令\(f'[n-j-i]=f[i+j]\)

再简单变形：\(f'[i]=f[n-i]\)

于是\(f'[]\)也构造出来了，这样两边都化为了卷积的形式，就可以用FFT求解了。

代码

#include<bits/stdc++.h>
#define M 400005
using namespace std;
int n;
struct Complex{
    double x,y;
    Complex(){}
    Complex(double _x,double _y):x(_x),y(_y){}
    Complex operator + (const Complex &res) const{
        return (Complex){x+res.x,y+res.y};  
    }
    Complex operator - (const Complex &res) const{
        return (Complex){x-res.x,y-res.y};  
    }
    Complex operator * (const Complex &res) const{
        return (Complex){x*res.x-y*res.y,x*res.y+y*res.x};  
    }
};
Complex A[M],B[M];
double pi=acos(-1.0);
void FFT(Complex *y,int n,int f){
    if(n==1)return;
    Complex L[n>>1],R[n>>1];
    for(int i=0;i<n;i+=2)L[i>>1]=y[i],R[i>>1]=y[i+1];
    FFT(L,n>>1,f);FFT(R,n>>1,f);
    Complex wn(cos(2*pi/n),f*sin(2*pi/n)),w(1,0);
    for(int i=0;i<(n>>1);i++,w=w*wn){
        y[i]=L[i]+w*R[i];
        y[i+(n>>1)]=L[i]-w*R[i];
    }
}
double q[M],b[M],C[M];
int main(){
    cin>>n;n--;
    for(int i=0;i<=n;i++)scanf("%lf",&q[i]);
    for(int i=0;i<=n;i++)A[i]=(Complex){q[i],0};
    for(int i=1;i<=n;i++){
        b[i]=1.0/i/i;
        B[i]=Complex(b[i],0);
    }
    int nn=n,mm=n;
    mm+=nn;
    for(nn=1;nn<=mm;nn<<=1);
    FFT(A,nn,1);FFT(B,nn,1);
    for(int i=0;i<=nn;i++)A[i]=A[i]*B[i];
    FFT(A,nn,-1);
    for(int i=0;i<=n;i++)C[i]=A[i].x/nn;
    for(int i=0;i<=nn;i++)A[i]=B[i]=Complex(0,0);
    for(int i=0;i<=n;i++)A[i]=Complex(q[n-i],0);
    for(int i=1;i<=n;i++)B[i]=Complex(b[i],0);
    FFT(A,nn,1);FFT(B,nn,1);
    for(int i=0;i<=nn;i++)A[i]=A[i]*B[i];
    FFT(A,nn,-1);
    for(int i=0;i<=n;i++)C[i]-=A[n-i].x/nn;
    for(int i=0;i<=n;i++)printf("%.3f\n",C[i]);
    return 0;
}