24 Sep 2016

BZOJ3527 力 FFT

题目大意:

给出$n$个数$q_i$,给出$F_i$的定义如下:

$F_i=\sum_{j<i} \frac {q_iq_j}{(i-j)^2}-\sum_{j>i} \frac{q_iq_j}{(i-j)^2}$

求所有$E_i=\frac {F_i} {q_i}$


题解:

怎么看这题都是非常的奇怪啊,在上面套了个$q_i$在求$E_i$的时候又除掉了。。

那么$E_i=\sum_{j<i} \frac {q_j}{(i-j)^2}-\sum_{j>i} \frac{q_j}{(i-j)^2}$

考虑构造卷积

揣摩一下$2$年前毛主力的心路历程:

这个卷积可以直接秒掉啊!

新开一个数组$p_j=\begin{cases} -(i-j)^{-2} && \text{j<n-i} \\ (i-j)^{-2} && \text{j>n-i} \\ 0 && \text{j=n-i}\end{cases}$

$E_i=\sum_{j=1}^n p_j q_{n-j}$

哎呀好像每次都要重新算怎么办

哎呀反正和$E_i$没什么关系只要得到那个卷积式子即可

为了让$q_i$不被实时计算,不如无脑加个$i$好了

于是变成了这样的东西:$E_{i}=\sum_{j=1}^{n}p_jq_{n+i-j}$

又当$j>n$时,$p_j=0$

所以式子还可以进一步扩展:

$E_i=\sum_{j=1}^{n+i}p_jq_{n+i-j}$

有点像卷积的形式了,把$E_i$换成$A_{i+n}$

于是就构成了卷积形式$A_{i+n}=\sum_{j=1}^{n+i}p_jq_{n+i-j}$

于是$FFT$一发可以得到$A$咯,再无脑减个$n$得到$E$

浙江OI真是简单


Code

#include<iostream>
#include<string.h>
#include<algorithm>
#include<stdio.h>
#include<math.h>
#include<complex>
using namespace std;
const int M=5e5+5;
typedef double db;
typedef complex<db> comp;
const db twpi=3.14159265358*2;
comp A[M],B[M],C[M];
int rev[M];
inline void DFT(comp *res,int n,int sg){
	for(int i=0;i<n;++i)
		if(rev[i]>i)swap(res[rev[i]+1],res[i+1]);
	for(int s=2;s<=n;s<<=1){
		comp wm(cos(twpi/s),sin(twpi/s)*sg);
		for(int i=1;i<=n;i+=s){
			comp w(1,0);
			for(int j=0;j<s>>1;++j,w*=wm){
				comp a=res[i+j],b=res[i+j+(s>>1)];
				res[i+j]=a+b*w,res[i+j+(s>>1)]=a-b*w;
			}
		}
	}
}
#define bug(x) cout<<#x<<"="<<x<<" "
#define debug(x) cout<<#x<<"="<<x<<endl
int main(){
	int n;cin>>n;
	db a;
	for(int i=1;i<=n;++i)
		scanf("%lf",&a),A[i]=a;
	int w=1,s=0;
	for(;w<=(n<<1);w<<=1,++s);
	for(int i=1;i<w;++i)
		rev[i]=rev[i>>1]>>1|(i&1)<<s-1;
	for(int i=1;i<=w;++i){
		if(i<n)B[i]=-1.0/((db)(n-i)*(n-i));
		else if(i>n)B[i]=1.0/((db)(n-i)*(n-i));
		else B[i]=0;
	}
	DFT(A,w,1),DFT(B,w,1);
	for(int i=1;i<=w;++i)
		C[i]=A[i]*B[i];
	DFT(C,w,-1);
	for(int i=1;i<=w;++i)C[i]/=w;
	for(int i=1;i<=n;++i)
		printf("%.3lf\n",C[i+n-1].real());
	return 0;
}

Tags:
Stats:
comments


Share:


Comments: