[BZOJ]3527: [Zjoi2014]力

Time Limit: 30 Sec  Memory Limit: 256 MB  Special Judge

Description

　　给出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}}$

　　令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

　　令$A(i)=q_{i}$，$B(i)=q_{n-i+1}$，$C(i)=1/i^{2}$，A与C卷积得到D，B与C卷积得到E，则答案$Ans(i)=D(i)-E(n-i+1)$，FFT卷积即可。

Code

#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
#define N 262144
struct cp
{
    double r,i;
    cp(double r=0,double i=0):r(r),i(i){}
    cp operator+(const cp&b){return cp(r+b.r,i+b.i);}
    cp operator-(const cp&b){return cp(r-b.r,i-b.i);}
    cp operator*(const cp&b){return cp(r*b.r-i*b.i,r*b.i+i*b.r);}
}a[N+5],b[N+5],c[N+5],w[2][N+5];
const double pi=acos(-1);
int r[N+5];
void init()
{
    cp g(cos(2*pi/N),sin(2*pi/N));int i,j,k;
    for(i=w[0][0].r=1;i<N;++i)w[0][i]=w[0][i-1]*g;
    for(i=w[1][0].r=1;i<N;++i)w[1][i]=w[0][N-i];
    for(i=j=0;i<N;r[++i]=j)for(k=N>>1;(j^=k)<k;k>>=1);
}
void fft(cp*x,int v)
{
    int i,j,k;
    for(i=1;i<N;++i)if(r[i]<i)swap(x[i],x[r[i]]);
    for(i=1;i<N;i<<=1)for(j=0;j<N;j+=i<<1)for(k=0;k<i;++k)
    {
        cp p=x[i+j+k]*w[v][N/(i<<1)*k];
        x[i+j+k]=x[j+k]-p;x[j+k]=x[j+k]+p;
    }
    if(v)for(i=0;i<N;++i)x[i].r/=N,x[i].i/=N;
}
int main()
{
    int n,i;
    scanf("%d",&n);
    for(i=1;i<=n;++i)scanf("%lf",&a[i].r),b[n-i+1]=a[i];
    for(i=1;i<=n;++i)c[i].r=1./i/i;
    init();fft(a,0);fft(b,0);fft(c,0);
    for(i=0;i<N;++i)a[i]=a[i]*c[i],b[i]=b[i]*c[i];
    fft(a,1);fft(b,1);
    for(i=1;i<=n;++i)printf("%.3lf
",a[i].r-b[n-i+1].r);
}