luogu P3338 [ZJOI2014]力

题目传送门：

https://www.luogu.org/problemnew/show/P3338

题意：

给定一个公式，求值。

思路：

考虑化简公式。

$F_j=\sum _{i<j}\frac{q_i\ast q_j}{(i-j{)}^2}-\sum _{i>j}\frac{q_i\ast q_j}{(i-j{)}^2}$

$F_j=\sum _{i=1}^{j-1}\frac{q_i\ast q_j}{(i-j{)}^2}-\sum _{i=j+1}^n\frac{q_i\ast q_j}{(i-j{)}^2}$

带入$E_j=\frac{F_j}{q_j}$，从而消除$q_j$，得：

$E_j=\sum _{i=1}^{j-1}\frac{q_i}{(i-j{)}^2}-\sum _{i=j+1}^n\frac{q_i}{(i-j{)}^2}$

令$A_i=\frac{1}{i^2}$，$q_i^{\prime }$为$q_i$反向取的序列，则有：

$E_j=\sum _{i=1}^{j-1}\frac{q_i}{(i-j{)}^2}-\sum _{i=1}^{j-1}\frac{q_i^{\prime }}{(i-j{)}^2}$

$E_j=\sum _{i=1}^{j-1}\frac{q_i}{(j-i{)}^2}-\sum _{i=1}^{j-1}\frac{q_i^{\prime }}{(j-i{)}^2}$

$E_j=\sum _{i=1}^{j-1}{q}_i\ast A_{j-i}-\sum _{i=1}^{j-1}\ast q_i^{\prime }\ast A_{j-i}$

容易发现这是一个类似于卷积的形式，便可以用$FFT$来加速。

代码：

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define LL long long
using namespace std;
const double pai=acos(-1.0);
struct comp
{
	double x,y;
	comp(double xx=0,double yy=0):x(xx),y(yy) {}
	friend comp operator+(const comp &x,const comp &y) {return comp(x.x+y.x,x.y+y.y);}
	friend comp operator-(const comp &x,const comp &y) {return comp(x.x-y.x,x.y-y.y);}
	friend comp operator*(const comp &a,const comp &b) {return comp(a.x*b.x-a.y*b.y,a.x*b.y+b.x*a.y);}
}a[500000],b[500000];
	int limit=1,n,l=0,ma=0;
	int r[500000];
	double d[500000][2],p[500000],f[500000];
void init()
{
	while(limit<=ma)
		limit<<=1,l++;
	for(int i=1;i<=limit;i++)
		r[i]=((r[i>>1]>>1)|((i&1)<<(l-1)));
}
void FFT(comp *now,int ty)
{
	for(int i=0;i<limit;i++)
		if(i<r[i]) swap(now[i],now[r[i]]);
	for(int mid=1;mid<limit;mid<<=1)
	{
		comp wn(cos(pai/mid),ty*sin(pai/mid));
		for(int j=0,R=(mid<<1);j<limit;j+=R)
		{
			comp w(1,0);
			for(int k=0;k<mid;k++,w=w*wn)
			{
				comp x=now[j+k],y=w*now[j+k+mid];
				now[j+k]=x+y;
				now[j+k+mid]=x-y;
			}
		}
	}
}
void work(int id)
{
	for(int i=0;i<=limit;i++)
		a[i]=(comp){0,0},b[i]=(comp){0,0};
	for(int i=1;i<=n;i++)
		a[i].x=d[i][id],b[i].x=p[i];
	FFT(a,1);
	FFT(b,1);
	for(int i=0;i<=limit;i++)
		a[i]=a[i]*b[i];
	FFT(a,-1);
	if(id==1)
		for(int i=0;i<=ma;i++)
			f[i]=a[i].x/limit;
	else
		for(int i=0;i<=ma;i++)
			f[i]-=a[n-i+1].x/limit;
}
int main()
{
	scanf("%d",&n);
	for(int i=1;i<=n;i++)
	{
		scanf("%lf",&d[i][1]);
		p[i]=1.0/i/i;
	}
	for(int i=1;i<=n;i++)
		d[i][2]=d[n-i+1][1];
	ma=(n<<1);
	init();
	work(1);
	work(2);
	for(int i=1;i<=n;i++)
		printf("%.3lf\n",f[i]);
}