BZOJ 3527: [Zjoi2014]力 [卷积][FFT]

最新推荐文章于 2019-08-26 14:30:00 发布

sillyf

最新推荐文章于 2019-08-26 14:30:00 发布

阅读量322

点赞数

分类专栏： FFT 文章标签： FFT

本文链接：https://blog.csdn.net/qq_35878547/article/details/79030756

版权

FFT 专栏收录该内容

2 篇文章 0 订阅

订阅专栏

题意

即求 $E_i=\sum_{j=1}^{i-1}(q_j/(j-i)^2)-\sum_{j=i+1}^{n}(q_j/(j-i)^2)$

题解

令 $f(x)=q_x\ ,g(x)=1/x^2$

则 $E_i=\sum_{j=1}^{i-1}f(j)g(j-i)-\sum_{j=i+1}^{n}f(j)g(j-i)$

发现这个式子很像卷积的形式,我们把 $q_i$ 下标存为 $[0,n-1]$ (以下记 $N=n-1$ )

然后令 $g(0)=0$ ,这样 $i==j$ 的情况也可以算进来了,并且对答案没有影响

分为两部分,左边那一块 $=\sum_{j=0}^{i}f(j)g(i-j)$ ,已经是卷积形式了

右边 $=\sum_{j=i}^{N}f(j)g(j-i)$ ,要转化首先让 $j$ 从零开始

右边 $=\sum_{j=0}^{N-i}f(j+i)g(j)$

构造 $F(x)=f(N-x)$

所以右边 $=\sum_{j=0}^{N-i}F(N-i-j)g(j)$

至此,两块都转为卷积形式,可以直接用fft加速计算

貌似考试时这题jz忘记给项数开两倍,导致爆炸了…

#include<cstdio>
#include<algorithm>
#include<cmath>
#define pi acos(-1.0)
#define N 2621450
using namespace std;

int n,m,R[N];
struct com{
    double r,i;
    com(){}
    com(double a,double b):r(a),i(b){}
}f[N],F[N],g[N];

inline com operator + (const com a,const com b){return com(a.r+b.r,a.i+b.i);}
inline com operator - (const com a,const com b){return com(a.r-b.r,a.i-b.i);}
inline com operator * (const com a,const com b){return com(a.r*b.r-a.i*b.i,a.i*b.r+a.r*b.i);}

inline void FFT(com *a,int f){
    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){
        com wn(cos(pi/k),f*sin(pi/k)),x,y;
        for(int i=0;i<n;i+=k<<1){
            com w(1,0);
            for(int j=0;j<k;++j,w=w*wn)
                x=a[i+j],y=w*a[i+j+k],a[i+j]=x+y,a[i+j+k]=x-y;
        }
    }
}

int main(){
//  freopen("3527.in","r",stdin);
//  freopen("3527.out","w",stdout);
    scanf("%d",&n);--n;
    for(int i=0;i<=n;++i){
        double x;
        scanf("%lf",&x);
        f[i].r=F[n-i].r=x;
        if(i) g[i].r=1.0/i/i;
    }
    m=n;int len=0;for(n=1;n<=(m<<1);n<<=1) ++len;
    for(int i=0;i<n;++i) R[i]=(R[i>>1]>>1)|((i&1)<<(len-1));
    FFT(f,1),FFT(F,1),FFT(g,1);
    for(int i=0;i<n;++i) f[i]=f[i]*g[i],F[i]=F[i]*g[i];
    FFT(f,-1),FFT(F,-1);
    for(int i=0;i<=m;++i) printf("%.3lf\n",(f[i].r-F[m-i].r)/n);
    return 0;
}