[BZOJ3527][Zjoi2014]力（FFT）

题目描述

传送门

题解

将式子化一下可以得出 E(i)=∑j=0iqj(i−j)2−∑j=inqj(i−j)2 只考虑第一个式子，令 f(n)=qn,g(n)=1n2 ，那么 F(i)=∑j=0if(i)g(i−j) 将序列 q 反转之后可以发现实际上第二个式子就是F(n−i) 正反求两遍卷积，然后答案相减即可

代码

#include<algorithm> #include<iostream> #include<cstring> #include<cstdio> #include<cmath> using namespace std; #define N 300005 const double pi=acos(-1.0); int n,L,R[N]; double q[N],c[N],d[N]; struct complex { double x,y; complex(double X=0,double Y=0) { x=X,y=Y; } }a[N],b[N]; complex operator + (complex a,complex b) {return complex(a.x+b.x,a.y+b.y);} complex operator - (complex a,complex b) {return complex(a.x-b.x,a.y-b.y);} complex operator * (complex a,complex b) {return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);} void FFT(complex a[N],int n,int opt) { for (int i=0;i<n;++i) if (i<R[i]) swap(a[i],a[R[i]]); for (int k=1;k<n;k<<=1) { complex wn=complex(cos(pi/k),opt*sin(pi/k)); for (int i=0;i<n;i+=(k<<1)) { complex w=complex(1,0); for (int j=0;j<k;++j,w=w*wn) { complex x=a[i+j],y=w*a[i+j+k]; a[i+j]=x+y,a[i+j+k]=x-y; } } } } void calc(int n,double *p) { int m=n+n;L=0;memset(R,0,sizeof(R)); for (n=1;n<=m;n<<=1) ++L; for (int i=0;i<n;++i) R[i]=(R[i>>1]>>1)|((i&1)<<(L-1)); FFT(a,n,1);FFT(b,n,1); for (int i=0;i<=n;++i) a[i]=a[i]*b[i]; FFT(a,n,-1); for (int i=0;i<=n;++i) p[i]=a[i].x/n; } int main() { scanf("%d",&n); for (int i=1;i<=n;++i) scanf("%lf",&q[i]); memset(a,0,sizeof(a));memset(b,0,sizeof(b)); a[0]=b[0]=complex(0,0); for (int i=0;i<=n;++i) a[i]=complex(q[i],0); for (int i=1;i<=n;++i) b[i]=complex(1.0/((double)i*(double)i),0); calc(n,c); memset(a,0,sizeof(a));memset(b,0,sizeof(b)); a[0]=b[0]=complex(0,0); for (int i=0;i<=n;++i) a[i]=complex(q[n-i],0); for (int i=1;i<=n;++i) b[i]=complex(1.0/((double)i*(double)i),0); calc(n,d); for (int i=1;i<=n;++i) printf("%.6lf\n",c[i]-d[n-i]); }