BZOJ 3527 [Zjoi2014]力(FFT)
题意:中文题
思路:我们对于这个公式拆成利两部分看,一部分是i<j,一部分是剩下的,我们在i<j的部分可以直接进行化简成卷积的形式,而对于大于的部分比较麻烦,相当于对所有的都替换成n-j,这样你两部分加起来的和还是定值,然后在做一次卷积即可
代码:(抄胡小兔爷的--->传送门)
#include <bits/stdc++.h> using namespace std; const int N = 1000005; const double PI=acos(-1.0); struct cp { double r,i; cp(){} cp(double _r,double _i) { r=_r;i=_i; } cp operator +(const cp a)const { return cp(a.r+r,a.i+i); } cp operator -(const cp a)const { return cp(r-a.r,i-a.i); } cp operator *(const cp a)const { return cp(r*a.r-i*a.i,r*a.i+i*a.r); } cp conj() { return cp(r,-i); } }; int n=1,m; double q[N],res[N]; cp f[N],f1[N],g[N],omg[N],inv[N]; void FFT_init(){ for(int i=0;i<n;i++){ omg[i]=cp(cos(2*PI*i/n),sin(2*PI*i/n)); inv[i]=omg[i].conj(); } } void fft(cp *a,cp *omg){ int lim=0; while((1 << lim) < n) lim++; for(int i=0;i<n;i++){ int t=0; for(int j=0;j<lim;j++) if(i >> j&1) t |= 1 << (lim - j - 1); if(i<t) swap(a[i],a[t]); } for(int l=2;l<=n;l*=2){ int m=l/2; for(cp *p=a;p!=a+n;p+=l){ for(int i=0;i<m;i++){ cp t=omg[n/l*i]*p[m+i]; p[m+i]=p[i]-t; p[i]=p[i]+t; } } } } int main(){ scanf("%d",&m); for(int i=0;i<m;i++){ scanf("%lf",&q[i]); f[i].r=q[i]; f1[m - 1 - i].r=q[i]; if(i) g[i].r=(1.0 / i / i); } while(n<2*m)n*=2; FFT_init(); fft(f,omg), fft(f1,omg), fft(g,omg); for(int i=0;i<n;i++) f[i]=f[i]*g[i],f1[i]=f1[i]*g[i]; fft(f,inv), fft(f1,inv); for(int i=0;i<m;i++){ res[i] = f[i].r/n - f1[m-1-i].r/n; } for(int i = 0; i < m; i++){ printf("%lf\n", res[i]); } return 0; }