题链:
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;
}