Description
给出n个数qi，给出Fj的定义如下：
令Ei=Fi/qi，求Ei.
Input
第一行一个整数n。
接下来n行每行输入一个数，第i行表示qi。
n≤100000，0<qi<1000000000
Output
n行，第i行输出Ei。与标准答案误差不超过1e-2即可。
Sample Input
5
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880
Sample Output
-16838672.693
3439.793
7509018.566
4595686.886
10903040.872
Solution
推式子
(displaystyle E_i=sum_{j<i}frac{q_j}{(i-j)^2}-sum_{j>i}frac{q_j}{(i-j)^2})
令 (g_i = egin{cases}&frac{1}{i^2}~~(i eq0)\&0 ~~~(i=0)end{cases})
那么
(displaystyle E_i=sum_{j=0}^{i-1}q_jg_{i-j}-sum_{j=i+1}^{n}q_jg_{i-j})
(displaystyle ~~~~~=sum_{j=0}^{i-1}q_jg_{i-j}-sum_{j=0}^{n-i-1}p_jg_{n-i-j}~~~~(p_i=g_{n-i}))
两个卷积形式，FFT就好了
#include<bits/stdc++.h>
#define ui unsigned int
#define ll long long
#define db double
#define ld long double
#define ull unsigned long long
const int MAXN=1<<19;
const db Pi=acos(-1.0);
int qn,n,m,cnt,rev[MAXN];
struct Complex{
db real,imag;
inline Complex operator + (const Complex &A) const {
return (Complex){real+A.real,imag+A.imag};
};
inline Complex operator - (const Complex &A) const {
return (Complex){real-A.real,imag-A.imag};
};
inline Complex operator * (const Complex &A) const {
return (Complex){real*A.real-imag*A.imag,imag*A.real+real*A.imag};
};
};
Complex q[MAXN],p[MAXN],g[MAXN];
template<typename T> inline void read(T &x)
{
T data=0,w=1;
char ch=0;
while(ch!='-'&&(ch<'0'||ch>'9'))ch=getchar();
if(ch=='-')w=-1,ch=getchar();
while(ch>='0'&&ch<='9')data=((T)data<<3)+((T)data<<1)+(ch^'0'),ch=getchar();
x=data*w;
}
template<typename T> inline void write(T x,char ch=' ')
{
if(x<0)putchar('-'),x=-x;
if(x>9)write(x/10);
putchar(x%10+'0');
if(ch!=' ')putchar(ch);
}
template<typename T> inline void chkmin(T &x,T y){x=(y<x?y:x);}
template<typename T> inline void chkmax(T &x,T y){x=(y>x?y:x);}
template<typename T> inline T min(T x,T y){return x<y?x:y;}
template<typename T> inline T max(T x,T y){return x>y?x:y;}
inline void FFT(Complex *A,int tp)
{
for(register int i=0;i<n;++i)
if(i<rev[i])std::swap(A[i],A[rev[i]]);
for(register int l=2;l<=n;l<<=1)
{
Complex wn=(Complex){cos(2*Pi/l),sin(tp*2*Pi/l)};
for(register int i=0;i<n;i+=l)
{
Complex w=(Complex){1,0};
for(register int j=0;j<(l>>1);++j)
{
Complex A1=A[i+j],A2=A[i+j+(l>>1)]*w;
A[i+j]=A1+A2;A[i+j+(l>>1)]=A1-A2;
w=w*wn;
}
}
}
}
int main()
{
read(qn);
m=qn+qn-1;
for(register int i=0;i<qn;++i)
{
scanf("%lf",&q[i].real);
p[qn-i].real=q[i].real;
if(!i)g[i].real=0.0;
else g[i].real=1.0/(1ll*i*i);
}
for(n=1;n<m;n<<=1)cnt++;
for(register int i=0;i<n;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(cnt-1));
FFT(q,1);FFT(p,1);FFT(g,1);
for(register int i=0;i<n;++i)q[i]=q[i]*g[i],p[i]=p[i]*g[i];
FFT(q,-1);FFT(p,-1);
for(register int i=0;i<qn;++i)printf("%f
",(q[i].real-p[qn-i].real)/n);
return 0;
}