题目地址:http://www.lydsy.com/JudgeOnline/problem.php?id=3527
题目大意:见原题。
算法讨论:
设A[i]=q[i],B[i]=1/(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即可。
Code:
/*
* Problem:3527
* Author:PYC
*/
#include <cstdio>
#include <cmath>
#define maxN 1000000
int n,N,L,dig[maxN+10],rev[maxN+10];
double X[maxN+10],ans[maxN+10];
const double pi=acos(-1);
struct complex{
double re,im;
complex(double r=0,double i=0){re=r,im=i;}
}x,y,a[maxN+10],b[maxN+10],c[maxN+10],A[maxN+10];
complex operator +(complex a,complex b){return complex(a.re+b.re,a.im+b.im);}
complex operator -(complex a,complex b){return complex(a.re-b.re,a.im-b.im);}
complex operator *(complex a,complex b){return complex(a.re*b.re-a.im*b.im,a.re*b.im+a.im*b.re);}
void FFT(complex a[],int f){
for (int i=0;i<N;++i) A[i]=a[rev[i]];
for (int i=0;i<N;++i) a[i]=A[i];
for (int i=2;i<=N;i<<=1){
complex wn(cos(2*pi/i),f*sin(2*pi/i));
for (int k=0;k<N;k+=i){
complex w(1);
for (int j=0;j<i/2;++j) x=a[k+j],y=w*a[k+j+i/2],a[k+j]=x+y,a[k+j+i/2]=x-y,w=w*wn;
}
}
if (f==-1) for (int i=0;i<N;++i) a[i].re/=N;
}
int main(){
scanf("%d",&n);
for (N=1;N<n;N<<=1,L++);N<<=1,L++;
for (int i=0;i<N;++i){
for (int t=i,len=0;t;t>>=1) dig[len++]=t&1;
for (int j=0;j<L;++j) rev[i]=rev[i]*2+dig[j];
}
for (int i=0;i<n;++i) scanf("%lf",&X[i]);
for (int i=0;i<n;++i) a[i]=complex(X[i]);
for (int i=1;i<n;++i) b[i]=complex(1./i/i);
FFT(a,1),FFT(b,1);
for (int i=0;i<N;++i) c[i]=a[i]*b[i];
FFT(c,-1);
for (int i=0;i<n;++i) ans[i]=c[i].re;
for (int i=0;i<N;++i) a[i]=b[i]=complex();
for (int i=0;i<n;++i) a[i]=complex(X[n-i-1]);
for (int i=1;i<n;++i) b[i]=complex(1./i/i);
FFT(a,1),FFT(b,1);
for (int i=0;i<N;++i) c[i]=a[i]*b[i];
FFT(c,-1);
for (int i=0;i<n;++i) ans[i]-=c[n-i-1].re;
for (int i=0;i<n;++i) printf("%0.9lf\n",ans[i]);
return 0;
}
By Charlie Pan
Mar 5,2014