BZOJ3527[Zjoi2014]力——FFT

题目描述

给出n个数qi，给出Fj的定义如下：

令Ei=Fi/qi，求Ei.

输入

第一行一个整数n。

接下来n行每行输入一个数，第i行表示qi。

n≤100000，0<qi<1000000000

输出

n行，第i行输出Ei。与标准答案误差不超过1e-2即可。

样例输入

5
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880

样例输出

-16838672.693
3439.793
7509018.566
4595686.886
10903040.872

将$q_{j}$移到等式左边，可以得到：$E_{k}=\sum\limits_{i<k}^{ }\frac{q_{i}}{(i-k)^2}-\sum\limits_{i>k}^{ }\frac{q_{i}}{(i-k)^2}$

将等式右边分成两部分看。

对于第一部分设$j=k-i$，可以得到：$E_{k}'=\sum\limits_{i+j==k}^{ }q_{i}*\frac{1}{j^2}$，设$F(i)=q_{i},G(i)=\frac{1}{i^2},H(i)=E_{i}'$，将$F$和$G$卷积即可得到$H$。

对于第二部分设$j=i-k$，可以得到：$E_{k}''=\sum\limits_{i-j==k}^{ }q_{i}*\frac{1}{j^2}$，即$E_{k}''=\sum\limits_{j+n-i+1==n-k+1}^{ }q_{i}*\frac{1}{j^2}$。

将$n-k+1$和$n-i+1$看作新的$k$和$i$，也就是将第一部分中的$H$和$F$翻转，然后再卷积即可。即$F(n-i+1)=q_{i},H(n-i+1)=E_{i}''$。

将两部分对应做减法即可。

#include<set>
#include<map>
#include<queue>
#include<stack>
#include<cmath>
#include<vector>
#include<cstdio>
#include<bitset>
#include<cstring>
#include<iostream>
#include<algorithm>
#define ll long long
using namespace std;
const double pi=acos(-1.0);
struct lty
{
	double x,y;
	lty(double X=0,double Y=0){x=X,y=Y;}
}f[400000],g[400000],h[400000];
int n,mask=1;
double ans[400000];
double q[400000];
lty operator +(lty a,lty b){return lty(a.x+b.x,a.y+b.y);}
lty operator -(lty a,lty b){return lty(a.x-b.x,a.y-b.y);}
lty operator *(lty a,lty b){return lty(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
void DFT(lty *a,int len)
{
	for(int i=0,k=0;i<len;i++)
	{
		if(i>k)
		{
			swap(a[i],a[k]);
		}
		for(int j=len>>1;(k^=j)<j;j>>=1);
	}
	for(int k=2;k<=len;k<<=1)
	{
		int t=k>>1;
		lty x(cos(pi/t),sin(pi/t));
		for(int i=0;i<len;i+=k)
		{
			lty w(1,0);
			for(int j=i;j<i+t;j++)
			{
				lty res=a[j+t]*w;
				a[j+t]=a[j]-res;
				a[j]=a[j]+res;
				w=w*x;
			}
		}
	}
}
void IDFT(lty *a,int len)
{
	for(int i=0,k=0;i<len;i++)
	{
		if(i>k)
		{
			swap(a[i],a[k]);
		}
		for(int j=len>>1;(k^=j)<j;j>>=1);
	}
	for(int k=2;k<=len;k<<=1)
	{
		int t=k>>1;
		lty x(cos(pi/t),-1.0*sin(pi/t));
		for(int i=0;i<len;i+=k)
		{
			lty w(1,0);
			for(int j=i;j<i+t;j++)
			{
				lty res=a[j+t]*w;
				a[j+t]=a[j]-res;
				a[j]=a[j]+res;
				w=w*x;
			}
		}
	}
	for(int i=0;i<len;i++)
	{
		a[i].x=a[i].x/len;
	}
}
int main()
{
	scanf("%d",&n);
	for(int i=1;i<=n;i++)
	{
		scanf("%lf",&q[i]);
	}
	while(mask<=(n<<1))
	{
		mask<<=1;
	}
	for(int i=1;i<=n;i++)
	{
		f[i].x=q[i];
	}
	for(int i=1;i<=n;i++)
	{
		g[i].x=(double)1/((double)i*i);
	}
	DFT(f,mask);
	DFT(g,mask);
	for(int i=0;i<mask;i++)
	{
		h[i]=g[i]*f[i];
	}
	IDFT(h,mask);
	for(int i=1;i<=n;i++)
	{
		ans[i]=h[i].x;
	}
	for(int i=0;i<mask;i++)
	{
		f[i].x=f[i].y=0;
		g[i].x=g[i].y=0;
	}
	for(int i=1;i<=n;i++)
	{
		f[i].x=q[n-i+1];
	}
	for(int i=1;i<=n;i++)
	{
		g[i].x=(double)1/((double)i*i);
	}
	DFT(f,mask);
	DFT(g,mask);
	for(int i=0;i<mask;i++)
	{
		h[i]=g[i]*f[i];
	}
	IDFT(h,mask);
	for(int i=1;i<=n;i++)
	{
		ans[i]-=h[n-i+1].x;
	}
	for(int i=1;i<=n;i++)
	{
		printf("%.3f\n",ans[i]);
	}
}