题链:
http://www.lydsy.com/JudgeOnline/problem.php?id=3527
题解:
FFT求卷积。
$$\begin{aligned}
E_i&=\frac{F_i}{q_i}\\
&=\sum_{k<i}\frac{q_k}{(i-k)^2}-\sum_{k>i}\frac{q_k}{(i-k)^2}\\
&=\sum_{k=1}^{n}{q_kP(i-k)}
\end{aligned}$$
其中,
$$P(x)=\begin{cases}
x^{-2}\;\;\quad x>0\\
0\;\;\quad\quad x=0\\
-x^{-2}\quad x<0\\
\end{cases}$$
如上,就构造出了一个卷积的形式,可以用FFT求出所有的$E_i$。
另外为了避免下标出现负数,所以E的下标和P的下标都加上n,(同时为了FFT的方便,把k从0开始枚举)即:
$$
E_i=E‘_{i+n}=\sum_{k=0}^{n-1}{q_kP‘(i+n-k)}\\
P‘(x)=\begin{cases}
(x-n)^{-2}\;\;\quad x>n\\
0\quad \qquad\qquad x=0\\
-(n-x)^{-2}\;\quad x<n\\
\end{cases}
$$
代码:
#include<bits/stdc++.h> #define MAXN 262145 using namespace std; const double Pi=acos(-1); struct Complex{ double real,image; Complex(double _real=0,double _image=0):real(_real),image(_image){} Complex operator - () const{return Complex(-real,-image);} friend Complex operator + (Complex A,Complex B){return Complex(A.real+B.real,A.image+B.image);} friend Complex operator - (Complex A,Complex B){return A+(-B);} friend Complex operator * (Complex A,Complex B){return Complex(A.real*B.real-A.image*B.image,A.image*B.real+A.real*B.image);} }; Complex Q[MAXN],P[MAXN]; int pos[MAXN]; void FFT(Complex *Y,int n,int sign){ for(int i=0;i<n;i++) if(i<pos[i]) swap(Y[i],Y[pos[i]]); for(int d=2;d<=n;d<<=1){ Complex dw(cos(2*Pi/d),sin(sign*2*Pi/d)),w,tmp; for(int i=0;w=Complex(1,0),i<n;i+=d) for(int k=i;k<i+d/2;w=w*dw,k++) tmp=w*Y[k+d/2],Y[k+d/2]=Y[k]-tmp,Y[k]=Y[k]+tmp; } } int main(){ int n,m,len; scanf("%d",&n); for(int i=0;i<n;i++) scanf("%lf",&Q[i].real); for(int k=1;k<n;k++) P[n-k].real=-1.0/(1.0*k*k),P[n+k].real=1.0/(1.0*k*k); m=2*n; for(n=1,len=0;n<m;n<<=1) len++; for(int i=0;i<n;i++) pos[i]=(pos[i>>1]>>1)|((i&1)<<(len-1)); FFT(Q,n,1); FFT(P,n,1); for(int i=0;i<n;i++) Q[i]=Q[i]*P[i]; FFT(Q,n,-1); for(int i=m/2;i<m;i++) printf("%.5lf\n",Q[i].real/n); return 0; }