欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页

[ZJOI2014]力 - FFT做卷积

程序员文章站 2022-05-22 19:19:45
...

Description
给出n个数qi,给出Fj的定义如下:
[ZJOI2014]力 - FFT做卷积
令Ei=Fi/qi,求Ei.
Input
第一行一个整数n。
接下来n行每行输入一个数,第i行表示qi。
n≤100000,0<qi<1000000000
Output
n行,第i行输出Ei。与标准答案误差不超过1e-2即可。

Sample Input
5

4006373.885184

15375036.435759

1717456.469144

8514941.004912

1410681.345880
Sample Output
-16838672.693

3439.793

7509018.566

4595686.886

10903040.872


看到这种卷积式子,我们应该想到FFT做卷积,然后我们式子就尽量往卷积上面化简。

Ej = sigma(i<j) qi / (i-j) / (i-j) - sigma(i>j) qi / (i-j) / (i-j)


我们首先要知道上面是卷积?

卷积的形式为: C[k] = sigma(i,0->k) ( A[i] * B[k-i] )

[ZJOI2014]力 - FFT做卷积[ZJOI2014]力 - FFT做卷积


然后直接FFT即可。

AC代码:

#pragma GCC optimize("-Ofast","-funroll-all-loops")
#include<bits/stdc++.h>
//#define int long long
using namespace std;
const int N=4e5+10;
const double PI=acos(-1.0);
int n,r[N],lim,l;
struct Complex{
	double x,y;
}a[N],b[N],c[N];
Complex operator + (Complex a,Complex b){return {a.x+b.x,a.y+b.y};}
Complex operator - (Complex a,Complex b){return {a.x-b.x,a.y-b.y};}
Complex operator * (Complex a,Complex b){return {a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
inline void FFT(Complex *a,int n,int k){
	for(int i=0;i<n;i++)	if(i<r[i])	swap(a[i],a[r[i]]);
	for(int mid=1;mid<n;mid<<=1){
		Complex wn={cos(2*PI/(mid<<1)),k*sin(2*PI/(mid<<1))};
		for(int i=0;i<n;i+=(mid<<1)){
			Complex w={1,0};
			for(int j=0;j<mid;j++,w=(w*wn)){
				Complex t0=a[i+j],t1=w*a[i+mid+j];
				a[i+j]=t0+t1;
				a[i+mid+j]=t0-t1;
			}
		}
	}
	if(k==-1)	for(int i=0;i<n;i++)	a[i].x/=n;
}
inline void init(){
	lim=1;	while(lim<=(n<<1))	lim<<=1,l++;
	for(int i=0;i<lim;i++)	r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
}
signed main(){
	cin>>n;
	for(int i=1;i<=n;i++){
		scanf("%lf",&a[i].x);
		c[n-i].x=a[i].x,b[i].x=(1.0/i/i);
	}
	init();
	FFT(a,lim,1),FFT(b,lim,1),FFT(c,lim,1);
	for(int i=0;i<lim;i++)	a[i]=a[i]*b[i],c[i]=c[i]*b[i];
	FFT(a,lim,-1),FFT(c,lim,-1);
	for(int i=1;i<=n;i++)	printf("%.3lf\n",a[i].x-c[n-i].x);
	return 0;
}
相关标签: FFT