【BZOJ】3527: [Zjoi2014]力 FFT

【参考】「ZJOI2014」力 - FFT by menci

【算法】FFT处理卷积

【题解】将式子代入后,化为Ej=Aj-Bj。

Aj=Σqi*[1/(i-j)^2],i=1~j-1。

令f(i)=qi,g(i)=1/i^2,定义f(0)=g(0)=0(方便卷积)。

Aj=Σf(i)*g(j-i),i=0~j-1,标准的卷积形式。

而对于Bj,将g反转后就是和为i+n-1的标准卷积形式了。

第一次FFT后,记得对a数组后半部分清零后再进行第二次FFT。

复杂度O(n log n)。

#include<cstdio>
#include<algorithm>
#include<complex>
#include<cmath>
using namespace std;
const int maxn=;
const double PI=acos(-);
complex<double>a1[maxn],a2[maxn];
int n;
double ans[maxn],b1[maxn],b2[maxn];
namespace fft{
complex<double>o[maxn],oi[maxn];
void init(int n){
for(int k=;k<n;k++)o[k]=complex<double>(cos(*PI*k/n),sin(*PI*k/n)),oi[k]=conj(o[k]);
}
void transform(complex<double>*a,int n,complex<double>*o){
int k=;
while((<<k)<n)k++;
for(int i=;i<n;i++){
int t=;
for(int j=;j<k;j++)if(i&(<<j))t|=(<<(k-j-));
if(i<t)swap(a[i],a[t]);
}
for(int l=;l<=n;l*=){
int m=l/;
for(complex<double>*p=a;p!=a+n;p+=l){
for(int i=;i<m;i++){
complex<double>t=o[n/l*i]*p[i+m];
p[i+m]=p[i]-t;
p[i]+=t;
}
}
}
}
void dft(complex<double>*a,int n){transform(a,n,o);}
void idft(complex<double>*a,int n){transform(a,n,oi);for(int i=;i<n;i++)a[i]/=n;}
}
void multply(int n){
int N=;
while(N<n+n)N*=;
for(int i=;i<N;i++)a1[i]=a2[i]=;
for(int i=;i<n;i++)a1[i].real(b1[i]),a2[i].real(b2[i]);
fft::init(N);fft::dft(a1,N);fft::dft(a2,N);
for(int i=;i<N;i++)a1[i]*=a2[i];
fft::idft(a1,N);
}
int main(){
scanf("%d",&n);n++;
b1[]=b2[]=;
for(int i=;i<n;i++){
scanf("%lf",&b1[i]);
b2[i]=1.0/i/i;
}
multply(n);
for(int i=;i<n;i++)ans[i]=a1[i].real();
for(int i=;i<n/;i++)swap(b2[i],b2[n-i-]);
multply(n);
for(int i=;i<n;i++){
ans[i]-=a1[i+n-].real();
printf("%.3lf\n",ans[i]);
}
return ;
}
上一篇:jQuery noConflict() 方法----与其他javaScript插件冲突时


下一篇:bzoj3527: [Zjoi2014]力 fft