1. 程式人生 > >BZOJ3527 || 洛谷P3338 [ZJOI2014]力【FFT】

BZOJ3527 || 洛谷P3338 [ZJOI2014]力【FFT】

Time Limit: 30 Sec
Memory Limit: 256 MB

Description

給出n個數qi,給出Fj的定義如下:
在這裡插入圖片描述
令Ei=Fi/qi,求Ei.

Input

第一行一個整數n。
接下來n行每行輸入一個數,第i行表示qi。
n≤100000,0<qi<1000000000

Output

n行,第i行輸出Ei。與標準答案誤差不超過1e-2即可。


題目分析

E i

= F i q i = j
&lt; i
q j ( i
j ) 2
j &gt; i q j ( j i ) 2 E_i=\frac{F_i}{q_i}=\sum_{j&lt;i}\frac{q_j}{(i-j)^2}-\sum_{j&gt;i}\frac{q_j}{(j-i)^2}
f ( i ) = q i , g ( i ) = 1 i 2 , h ( i ) = q n i + 1 f(i)=q_i,g(i)=\frac{1}{i^2},h(i)=q_{n-i+1}

則有 E i = j = 1 i 1 f ( j ) g ( i j ) j = i + 1 n h ( n j + 1 ) g ( j i ) E_i=\sum_{j=1}^{i-1}f(j)*g(i-j)-\sum_{j=i+1}^nh(n-j+1)*g(j-i)
發現這個式子兩部分都是卷積形式,故FFT優化即可


#include<iostream>
#include<cmath>
#include<algorithm>
#include<map>
#include<cstring>
#include<cstdio>
using namespace std;
typedef long long lt;
typedef double dd;
 
int read()
{
    int f=1,x=0;
    char ss=getchar();
    while(ss<'0'||ss>'9'){if(ss=='-')f=-1;ss=getchar();}
    while(ss>='0'&&ss<='9'){x=x*10+ss-'0';ss=getchar();}
    return f*x;
}

const dd Pi=acos(-1.0);
const int maxn=400010;
int n,m;
struct complex{
    dd x,y;
    complex(dd _x=0,dd _y=0){ x=_x; y=_y;}
}A[maxn],B[maxn],C[maxn];
int lim=1,L,R[maxn];

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

void FFT(complex* a,int opt)
{
    for(int i=0;i<lim;++i)
    if(i<R[i]) swap(a[i],a[R[i]]);
    
    for(int i=1;i<lim;i<<=1)
    {
        complex wn(cos(Pi/i),opt*sin(Pi/i));
        for(int j=0;j<lim;j+=(i<<1))
        {
            complex w(1,0);
            for(int k=0;k<i;++k)
            {
                complex nx=a[j+k],ny=w*a[i+j+k];
                a[j+k]=nx+ny;
                a[i+j+k]=nx-ny;
                w=w*wn;
            }
        }
    }
}

int main()
{
    n=read();
    for(int i=1;i<=n;++i) 
    {
        scanf("%lf",&A[i].x);
        B[n-i+1].x=A[i].x;
		C[i].x=1.0/(dd)i/(dd)i;
    } 
    
    while(lim<=(n<<1)) lim<<=1,L++;
    for(int i=0;i<lim;++i)
    R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
    
    FFT(A,1); FFT(B,1); FFT(C,1);
    for(int i=0;i<=lim;++i) 
    A[i]=A[i]*C[i],B[i]=B[i]*C[i];
    
    FFT(A,-1); FFT(B,-1);
    for(int i=1;i<=n;++i)
    printf("%.3lf\n",(A[i].x-B[n-i+1].x)/(dd)lim);
    return 0;
    
}