P3338 [ZJOI2014]力

题意

首先化简式子:
(E_j=frac{F_j}{q_j})
(=frac{sumlimits_{i=1}^{j-1}frac{q_i*q_j}{(i-j)^2}-sumlimits_{i=j+1}^{n}frac{q_i*q_j}{(i-j)^2}}{q_j})
(=sumlimits_{i=1}^{j-1}frac{q_i}{(i-j)^2}-sumlimits_{i=j+1}^{n}frac{q_i}{(i-j)^2})
(f_i=q_i,g_i=frac{1}{i^2})
(=sumlimits_{i=1}^{j-1}f_i*g_{j-i}-sumlimits_{i=j+1}^nf_i*g_{j-i})
显然前面那个式子就是个卷积的形式,可以直接(FFT),考虑后面那个式子:
(sumlimits_{i=j+1}^nf_i*g_{j-i})
(f'(x))表示(f(x))翻转后的多项式:
(=sumlimits_{i=1}^{j-1}f'_i*g_{j-i})

于是分别对两式做(FFT),最后将第二个翻转后相减即可。

code:

#include<bits/stdc++.h>
using namespace std;
const int maxn=1e6+10;
const double Pi=acos(-1.0);
int n,lim=1,len;
int pos[maxn];
double q[maxn];
struct cplx{double x,y;}f[maxn],g[maxn],h[maxn],a[maxn],b[maxn];
cplx operator+(cplx a,cplx b){return (cplx){a.x+b.x,a.y+b.y};}
cplx operator-(cplx a,cplx b){return (cplx){a.x-b.x,a.y-b.y};}
cplx operator*(cplx a,cplx b){return (cplx){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
inline void FFT(cplx* a,int op)
{
	for(int i=0;i<lim;i++)if(i<pos[i])swap(a[i],a[pos[i]]);
	for(int l=1;l<lim;l<<=1)
	{
		cplx wn=(cplx){cos(Pi/l),op*sin(Pi/l)};
		for(int i=0;i<lim;i+=l<<1)
		{
			cplx w=(cplx){1,0};
			for(int j=0;j<l;j++,w=w*wn)
			{
				cplx x=a[i+j],y=w*a[i+l+j];
				a[i+j]=x+y,a[i+l+j]=x-y;
			}
		}
	}
	if(op==1)return;
	for(int i=0;i<lim;i++)a[i].x/=lim;
}
int main()
{
	scanf("%d",&n);
	for(int i=1;i<=n;i++)scanf("%lf",&q[i]);
	for(int i=1;i<=n;i++)f[i].x=h[i].x=q[i],g[i].x=1.0/(1.0*i*i);
	reverse(h+1,h+n+1);
	while(lim<(n<<1))lim<<=1,len++;
	for(int i=0;i<lim;i++)pos[i]=(pos[i>>1]>>1)|((i&1)<<(len-1));
	FFT(f,1);FFT(g,1);FFT(h,1);
	for(int i=0;i<lim;i++)a[i]=f[i]*g[i],b[i]=h[i]*g[i];
	FFT(a,-1);FFT(b,-1);
    reverse(b+1,b+n+1);
	for(int i=1;i<=n;i++)printf("%.3lf
",a[i].x-b[i].x);
	return 0;
}
原文地址:https://www.cnblogs.com/nofind/p/12133211.html