BZOJ 3527: [Zjoi2014]力

题目地址：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