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\%\)

對於\(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