【BZOJ】3527: [Zjoi2014]力 FFT

【参考】「ZJOI2014」力 - FFT by menci

【算法】FFT处理卷积

【题解】将式子代入后，化为Ej=Aj-Bj。

Aj=Σqi*[1/(i-j)^2]，i=1~j-1。

令f(i)=qi，g(i)=1/i^2，定义f(0)=g(0)=0（方便卷积）。

Aj=Σf(i)*g(j-i)，i=0~j-1，标准的卷积形式。

而对于Bj，将g反转后就是和为i+n-1的标准卷积形式了。

第一次FFT后，记得对a数组后半部分清零后再进行第二次FFT。

复杂度O(n log n)。

#include<cstdio>
#include<algorithm>
#include<complex>
#include<cmath>
using namespace std;
const int maxn=300010;
const double PI=acos(-1);
complex<double>a1[maxn],a2[maxn];
int n;
double ans[maxn],b1[maxn],b2[maxn];
namespace fft{
    complex<double>o[maxn],oi[maxn];
    void init(int n){
        for(int k=0;k<n;k++)o[k]=complex<double>(cos(2*PI*k/n),sin(2*PI*k/n)),oi[k]=conj(o[k]);
    }
    void transform(complex<double>*a,int n,complex<double>*o){
        int k=0;
        while((1<<k)<n)k++;
        for(int i=0;i<n;i++){
            int t=0;
            for(int j=0;j<k;j++)if(i&(1<<j))t|=(1<<(k-j-1));
            if(i<t)swap(a[i],a[t]);
        }
        for(int l=2;l<=n;l*=2){
            int m=l/2;
            for(complex<double>*p=a;p!=a+n;p+=l){
                for(int i=0;i<m;i++){
                    complex<double>t=o[n/l*i]*p[i+m];
                    p[i+m]=p[i]-t;
                    p[i]+=t;
                }
            }
        }
    }
    void dft(complex<double>*a,int n){transform(a,n,o);}
    void idft(complex<double>*a,int n){transform(a,n,oi);for(int i=0;i<n;i++)a[i]/=n;}
}
void multply(int n){
    int N=1;
    while(N<n+n)N*=2;
    for(int i=0;i<N;i++)a1[i]=a2[i]=0;
    for(int i=0;i<n;i++)a1[i].real(b1[i]),a2[i].real(b2[i]);
    fft::init(N);fft::dft(a1,N);fft::dft(a2,N);
    for(int i=0;i<N;i++)a1[i]*=a2[i];
    fft::idft(a1,N);
}
int main(){
    scanf("%d",&n);n++;
    b1[0]=b2[0]=0;
    for(int i=1;i<n;i++){
        scanf("%lf",&b1[i]);
        b2[i]=1.0/i/i;
    }
    multply(n);
    for(int i=1;i<n;i++)ans[i]=a1[i].real();
    for(int i=0;i<n/2;i++)swap(b2[i],b2[n-i-1]);
    multply(n);
    for(int i=1;i<n;i++){
        ans[i]-=a1[i+n-1].real();
        printf("%.3lf\n",ans[i]);
    }
    return 0;
}