[ZJOI2014][JZOJ3617]力

题目大意

给出$n$个数$q_i$，设

$$F_j=\sum_{i<j}\frac{q_iq_j}{(i-j)^2}-\sum_{i>j}\frac{q_iq_j}{(i-j)^2}$$
令 $E_i=\frac{F_i}{q_i}$试求 $E_i$。
$1\le n\le100000$， $0\lt q_i\lt1000000000$

---

题目分析

这是我第一道$\mathrm{FFT}$，搞了很久。
我们简化一下

$$\begin{align} E_j&=\frac{\sum_{i<j}\frac{q_iq_j}{(i-j)^2}-\sum_{i>j}\frac{q_iq_j}{(i-j)^2}}{q_j}\\ &=\sum_{i<j}\frac{q_i}{(i-j)^2}-\sum_{i>j}\frac{q_i}{(i-j)^2} \end{align}$$
显然可以分成前后两部分算，第二部分直接将 $q$数组倒过来就和第一部分一样了。
那我们怎么做第一部分呢？令数组$b_i=\frac{1}{i^2}(i\not=0)$， $b_0=0$，那么我们就可以将上式变为标准点积
$$Ans_k=\sum_{i=0}^{k}q_ib_{k-i}$$
然后上 $\mathrm{FFT}$就行了。
时间复杂度 $\mathrm O(nlog_2^n)$。

---

代码实现

#include <iostream>
#include <cstring>
#include <cstdio>
#include <cmath>

using namespace std;

const int N=262144;
const double pi=acos(-1);

struct C
{
    double x,y;
    C (double x0=0,double y0=0){x=x0,y=y0;}
};

C operator+(C a,C b){return C(a.x+b.x,a.y+b.y);}
C operator-(C a,C b){return C(a.x-b.x,a.y-b.y);}
C operator*(C a,C b){return C(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
C operator*(C a,double b){return C(a.x*b,a.y*b);}
C operator/(C a,C b){return C((a.x*b.x+a.y*b.y)/(b.x*b.x-b.y*b.y),(a.y*b.x-a.x*b.y)/(b.x*b.x-b.y*b.y));}
C operator/(C a,double b){return C(a.x/b,a.y/b);}

C f[N],g[N],f_[N],r1[N],r2[N],t[N];
double q[N],ans[N];
int id[N];
int n,n0,l;

void preparation()
{
    l=(int)ceil(log(n0*2)/log(2));
    n=1<<l;
    for (int i=0;i<n;i++)
        for (int j=i,k=0;k<l;j/=2,k++)
            id[i]=(id[i]<<1)+(j&1);
}

void DFT(C *a,int sig)
{
    for (int i=0;i<n;i++)
        t[id[i]]=a[i];
    for (int m=2;m<=n;m*=2)
    {
        int h=m>>1;
        for (int i=0;i<h;i++)
        {
            C w(cos(i*sig*pi/h),sin(i*sig*pi/h));
            for (int j=i;j<n;j+=m)
            {
                C u=t[j],v=w*t[j+h];
                t[j]=u+v;
                t[j+h]=u-v;
            }
        }
    }
    for (int i=0;i<n;i++)
        a[i]=t[i];
    if (sig==-1)
        for (int i=0;i<n;i++)
            a[i]=a[i]/n;
}

void FFT(C *a,C *b,C *c)
{
    DFT(a,1);
    DFT(b,1);
    for (int i=0;i<n;i++)
        c[i]=a[i]*b[i];
    DFT(c,-1);
}

void solve()
{
    g[0]=C(0,0);
    for (int i=1;i<n0;i++)
        g[i]=C(1.0/i/i,0);
    FFT(f,g,r1);
    memset(g,0,sizeof g);
    for (int i=1;i<n0;i++)
        g[i]=C(1.0/i/i,0);
    FFT(f_,g,r2);
    for (int i=0;i<n0;i++)
        ans[i]=r1[i].x-r2[n0-i-1].x;
}

int main()
{
    freopen("force.in","r",stdin);
    freopen("force.out","w",stdout);
    scanf("%d",&n0);
    for (int i=0;i<n0;i++)
    {
        scanf("%lf",&q[i]);
        f[i]=C(q[i],0);
        f_[n0-i-1]=C(q[i],0);
    }
    preparation();
    solve();
    for (int i=0;i<n0;i++)
        printf("%.4lf\n",ans[i]);
    fclose(stdin);
    fclose(stdout);
    return 0;
}