P3338 [ZJOI2014]力

题目描述

给出n个数qi，给出Fj的定义如下：

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

令\(E_i=\frac{F_i}{q_i}\)，求\(E_i\).

输入输出格式

输入格式：

第一行一个整数\(n\)。

接下来\(n\)行每行输入一个数，第\(i\)行表示\(q_i\)。

输出格式：

\(n\)行，第\(i\)行输出\(E_i\)。

与标准答案误差不超过\(10^{-2}\)即可。

说明

对于\(30\%\)的数据，\(n\le1000\)。

对于\(50\%\)的数据，\(n\le60000\)。

对于\(100\%\)的数据，\(n\le100000\)，\(0<qi<1000000000\)。

[spj 0.01]

---

鉴于本*一开始啥也看不出来还是写一写好了。

把两边分开考虑。令\(f(x)=\frac{1}{x^2},g(x)=q_x\)，然后发现就是\(FFT\)的形式，看成系数做就可以了。

右边转过来以后同理。

有一些细节比如

• 1/(i*i)会被卡精度，要写1/i/i
• \(f(0)=g(0)=0\)

---

Code:

#include <cstdio>

#include <cmath>

#include <algorithm>

const int N=(1<<18)+10;

struct complex

{

    double x,y;

    complex(){}

    complex(double x,double y){this->x=x,this->y=y;}

    complex friend operator +(complex n1,complex n2){return complex(n1.x+n2.x,n1.y+n2.y);}

    complex friend operator -(complex n1,complex n2){return complex(n1.x-n2.x,n1.y-n2.y);}

    complex friend operator *(complex n1,complex n2){return complex(n1.x*n2.x-n1.y*n2.y,n1.x*n2.y+n1.y*n2.x);}

}a[N],b[N],tmpx,tmpy,wn,w;

const double pi=3.1415926535897632;

int n,turn[N],len=1,L=-1;

double out[N],q[N];

void FFT(complex *a,int typ)

{

    for(int i=0;i<len;i++)

        if(i<turn[i])

            std::swap(a[i],a[turn[i]]);

    for(int le=1;le<len;le<<=1)

    {

        wn=complex(cos(pi/le),typ*sin(pi/le));

        for(int p=0;p<len;p+=le<<1)

        {

            w=complex(1,0);

            for(int i=p;i<p+le;i++,w=w*wn)

            {

                tmpx=a[i],tmpy=w*a[i+le];

                a[i]=tmpx+tmpy,a[i+le]=tmpx-tmpy;

            }

        }

    }

}

int main()

{

    scanf("%d",&n);

    for(int i=1;i<=n;i++) scanf("%lf",q+i);

    for(int i=1;i<=n;i++) a[i].x=q[i],b[i].x=1.0/i/i;

    while(len<=n+1<<1) len<<=1,++L;

    for(int i=0;i<len;i++) turn[i]=turn[i>>1]>>1|(i&1)<<L;

    FFT(a,1),FFT(b,1);

    for(int i=0;i<len;i++) a[i]=a[i]*b[i];

    FFT(a,-1);

    for(int i=1;i<=n;i++) out[i]+=a[i].x/len;

    for(int i=0;i<len;i++) a[i]=complex(0,0);

    for(int i=1;i<=n;i++)

        a[n+1-i]=complex(q[i],0);

    FFT(a,1);

    for(int i=0;i<len;i++) a[i]=a[i]*b[i];

    FFT(a,-1);

    for(int i=1;i<=n;i++) out[n+1-i]-=a[i].x/len;

    for(int i=1;i<=n;i++) printf("%lf\n",out[i]);

    return 0;

}

---

2018.12.3