JZOJ 3617. 【ZJOI2014】力

Description

Input

输入文件包含一个整数n，接下来n行每行输入一个数，第i行表示qi。

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

Data Constraint

对于30%的数据，n≤1000。

对于50%的数据，n≤60000。

对于100%的数据，n≤100000，0

Solution

• 经典FFT……

• 先推式子，化简一下 $E_i$ ，得：

  $$E_j=\sum_{i<j}\frac{q_i}{(j-i)^2}-\sum_{i>j}\frac{q_i}{(i-j)^2}$$

• 可以发现前面和后面是差不多的，我们可以分开算。

• 现在考虑前面，有：

  $$P_j=\sum_{i<j}\frac{q_i}{(j-i)^2}$$

• 令 $r_i=\frac{1}{i^2}$ ，则有标准数论卷积的形式：

  $$P_j=\sum_{i<j}q_i*r_{j-i}$$

• 如此一来就可以直接用FFT加速计算了，时间复杂度 $O(N\ log\ N)$ 。

• 而后面部分就相当于把 q <script type="math/tex" id="MathJax-Element-533">q</script> 数组反过来做一遍即可。

Code

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cctype>
#include<cmath>
using namespace std;
const int N=1e5+5;
const double Pi=acos(-1.0);
struct comp
{
    double r,i;
    comp(){}
    comp(double rr,double ii){r=rr,i=ii;}
    comp operator +(const comp &x)const
    {
        return comp(r+x.r,i+x.i);
    }
    comp operator -(const comp &x)const
    {
        return comp(r-x.r,i-x.i);
    }
    comp operator *(const comp &x)const
    {
        return comp(r*x.r-i*x.i,i*x.r+r*x.i);
    }
}f1[N<<2],f2[N<<2],f3[N<<2];
int n,m;
int rev[N<<2];
double q[N];
inline double read()
{
    double X=0,Y=1.0; int w=0; char ch=0;
    while(!isdigit(ch)) w|=ch=='-',ch=getchar();
    while(isdigit(ch)) X=X*10+(ch^48),ch=getchar();
    ch=getchar();
    while(isdigit(ch)) X+=(Y/=10)*(ch^48),ch=getchar();
    return w?-X:X;
}
inline void FFT(comp *y,int ff)
{
    for(int i=0;i<m;i++)
        if(i<rev[i]) swap(y[i],y[rev[i]]);
    for(int h=2;h<=m;h<<=1)
    {
        comp wn(cos(2*Pi/h),ff*sin(2*Pi/h));
        for(int i=0;i<m;i+=h)
        {
            comp w(1,0);
            for(int k=i;k<i+h/2;k++)
            {
                comp u=y[k],t=w*y[k+h/2];
                y[k]=u+t;
                y[k+h/2]=u-t;
                w=w*wn;
            }
        }
    }
    if(ff==-1) for(int i=0;i<m;i++) y[i].r/=m;
}
int main()
{
    scanf("%d",&n);
    for(int i=1;i<=n;i++) q[i]=read();
    int l=0;
    for(m=1;m<=n<<1;m<<=1) l++;
    for(int i=0;i<m;i++) rev[i]=rev[i>>1]>>1|(i&1)<<l-1;
    for(int i=1;i<=n;i++) f1[i]=comp(q[i],0);
    for(int i=1;i<=n;i++) f2[i]=comp(q[n-i+1],0);
    for(int i=1;i<=n;i++) f3[i]=comp(1.0/i/i,0);
    FFT(f1,1),FFT(f2,1),FFT(f3,1);
    for(int i=0;i<m;i++) f1[i]=f1[i]*f3[i];
    for(int i=0;i<m;i++) f2[i]=f2[i]*f3[i];
    FFT(f1,-1),FFT(f2,-1);
    for(int i=1;i<=n;i++) printf("%.3lf\n",f1[i].r-f2[n-i+1].r);
    return 0;
}