此文中介绍了FFT的基本套路:
本题便是对FFT基本套路的又一次实践,即翻转一个串,把对应相乘被转化成卷积。
题目描述
给出 n 个数
qi , Fj 的定义如下:
Fj=∑i<jqiqj(i−j)2−∑i>jqiqj(i−j)2令: Ei=Fi/qi ,求 Ei 。
输入格式
第一行一个数n,接下来n行每行一个数,第i行的数为 qi 。
输出格式
n行,第i行的数表示 Ei 。
样例输入
5
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880
样例输出
-16838672.693
3439.793
7509018.566
4595686.886
10903040.872
解题思路
把 Fj 中的所有 qj 都提出来,有:
所以有:
把这个式子拆成两半,令:
则有:
先看 f(i) : f(j)=∑j−1i=0qi(i−j)2 ,不妨令 i=j−i (调整循环顺序,结果不变),有:
定义向量A,B,使得: Ai=qi , Bi=1i2(i∈1..n) ,特别地, B0=0 。则有:
令向量C为,向量A、B的卷积。即: C=A×B ,有:
由此证明:
然后考虑 g(j) : g(j)=∑ni=j+1qi(i−j)2 不妨令 i=j−i (同上文,调整循环顺序,结果不变),有:
i小于零看着太别扭,再令 i=−i ,有:
此时,我们发现,i与j+i的差是一个定值。
定义一个向量 B′ ,使得 B′i=Bn−i ,即对B向量进行翻转。那么就有:
这次,不难看出,n-i与j+i的和是一个定值了。
定义向量 C′ ,有 C′=A×B′ ,即:
那么:
令 i=j+i ,有:
当 i<0 时, n−i>n , B′n−i=0 ,可以不考虑;当 i>n−j 时, j+i>n , Aj+i=0 ,也可以不考虑,所以:
结合上文中所述,我们可以利用FFT以 O(nlgn) 的时间复杂度计算卷积 C=A×B 和 C′=A×B′ 。输出的结果就是 Ei=Ci−C′n+i 。
代码
#include<algorithm>
#include<cmath>
#include<complex>
using namespace std;
typedef complex<double>cd;
#include<stdio.h>
#include<stdlib.h>
#define maxl ((100000+3)*4)
cd A[maxl],B[maxl],invB[maxl];
cd f[maxl],g[maxl];
double ans[maxl];
int rev[maxl];
void get_rev(int bit){
for(int i=0;i<(1<<bit);i++){
rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
}
}
void fft(cd* a,int n,int dft){
for(int i=0;i<n;i++){
if(i<rev[i])
swap(a[i],a[rev[i]]);
}
for(int step=1;step<n;step<<=1){
cd wn=exp(cd(0,dft*M_PI/step));
for(int j=0;j<n;j+=step<<1){
cd wnk(1,0);
for(int k=j;k<j+step;k++){
cd x=a[k];
cd y=wnk*a[k+step];
a[k]=x+y;
a[k+step]=x-y;
wnk*=wn;
}
}
}
if(dft==-1){
for(int i=0;i<n;i++)
a[i]/=n;
}
}
int main(){
int n;
scanf("%d",&n);
for(int i=0;i<n;i++){
double x;
scanf("%lf",&x);
A[i]=x;
}
for(int i=1;i<=n;i++)B[i]=1.0/i/i;
for(int i=0;i<=n;i++)invB[i]=B[n-i];
int bit,s;
for(bit=1,s=2;(1<<bit)<(n+(n+1))-1;bit++)s<<=1;
get_rev(bit);
fft(A,s,1);fft(B,s,1);fft(invB,s,1);
for(int i=0;i<s;i++){
f[i]=A[i]*B[i];
//printf("f[%d]=(%lf,%lf)\n",i,f[i].real(),f[i].imag());
g[i]=A[i]*invB[i];
}
fft(f,s,-1);fft(g,s,-1);
for(int i=0;i<n;i++){
ans[i]=f[i].real()-g[n+i].real();
printf("%.3lf\n",ans[i]);
}
return 0;
}