P3338 [ZJOI2014]力

蒟蒻第一次做不是板子(?不是吗)的fft题(竟然很流畅地50minA了)

首先化简式子

$E_{i}=\sum_{j<i}^{ }\frac{q_{j}}{(i-j)^{2}}-\sum_{j>i}^{ }\frac{q_{j}}{(i-j)^{2}}$

然后可以设数列{$q_{i}$}与{$a_{i}$}

{$a_{i}$}的定义可以看代码,不太好描述

然后$ans_{i}=\sum_{k+j=n+i}q_{j}*a_{k}$(不知道写法有没有问题...)

然后可以发现这是卷积的形式

然后fft就可以快乐卷积辣

code:

#include<bits/stdc++.h>
using namespace std;
const int N=1000000;
const double pi=acos(-1);
struct comp{
    double x,y;
    comp(double xx=0,double yy=0){
        x=xx,y=yy;
    }
};
comp operator +(comp a,comp b){
    return comp(a.x+b.x,a.y+b.y);
}
comp operator -(comp a,comp b){
    return comp(a.x-b.x,a.y-b.y);
}
comp operator *(comp a,comp b){
    return comp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
}
comp a[N],b[N];
int n,L,lim=1,r[N];
void fft(comp *f,int op){
    for(int i=0;i<lim;i++){
        if(i<r[i])swap(f[i],f[r[i]]);
    }
    for(int p=2;p<=lim;p<<=1){
        int len=p>>1;
        comp tem(cos(pi/len),op*sin(pi/len));
        for(int k=0;k<lim;k+=p){
            comp buf(1,0);
            for(int j=k;j<k+len;j++){
                comp tt=buf*f[j+len];
                f[j+len]=f[j]-tt;
                f[j]=f[j]+tt;
                buf=buf*tem;
            }
        }
    }
}
int main(){
    cin>>n;
    for(int i=0;i<n;i++){
        double x;
        scanf("%lf",&x);
        b[i].x=x;
    }
    for(int i=-n+1;i<=n-1;i++){
        if(i<0)
            a[i+n-1].x=-1.0/i/i;
        else if(i==0)
            a[i+n-1].x=0;
        else
            a[i+n-1].x=1.0/i/i;
    }
    while((1<<L)<(3*n-3))L++,lim<<=1;
    for(int i=0;i<lim;i++){
        r[i]=(r[i>>1]>>1)|((i&1)<<(L-1));
    }
    fft(a,1);
    fft(b,1);
    for(int i=0;i<lim;i++)a[i]=a[i]*b[i];
    fft(a,-1);
    for(int i=n-1;i<=2*n-2;i++){
        printf("%.3f\n",a[i].x/lim);
    }
    return 0;
}

 

posted @ 2019-07-19 14:58  天才美少女雪乃  阅读(160)  评论(0编辑  收藏  举报