1. 程式人生 > >[Luogu]P3338 [ZJOI2014]力(FFT)

[Luogu]P3338 [ZJOI2014]力(FFT)

題目描述

給出\(n\)個數\(q_i\),給出\(F_j\)的定義如下:

\(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=F_i/q_i\),求\(E_i\).

輸入輸出格式

輸入格式:

第一行一個整數n。

接下來n行每行輸入一個數,第i行表示qi。

輸出格式:

n行,第i行輸出Ei。

與標準答案誤差不超過1e-2即可。

輸入輸出樣例

輸入樣例#1: 複製

5
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880

輸出樣例#1: 複製

-16838672.693
3439.793
7509018.566
4595686.886
10903040.872

說明

對於30%的資料,n≤1000。

對於50%的資料,n≤60000。

對於100%的資料,n≤100000,0<qi<1000000000。

[spj 0.01]

題解

一道對於剛學FFT的人不錯的題目。
完全可以自己手推。
搞了我晚自習半個小時才推出來...作業都沒寫完

對於這個式子
\(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
所以我們首先消掉qi
所以變成了這樣
\(E_j = \sum_{i<j}\frac{q_i}{(i-j)^2 }-\sum_{i>j}\frac{q_i}{(i-j)^2 }\)


好,接下來我們發現暴力O(n^2)就可以求出來了。
怎麼轉換到FFT裡面去?
先O(n)處理出(i-j)^2.
這個時候我們列出樣例。
對於E1=

4006373.885184*(1-1)^2 +
15375036.435759*(1-2)^2 -
1717456.469144*(1-3)^2 -
8514941.004912*(1-4)^2 -
1410681.345880*(1-5)^2 -

然後加起來
對於E2=

4006373.885184*(2-1)^2 +
15375036.435759*(2-2)^2 +
1717456.469144*(2-3)^2 -
8514941.004912*(2-4)^2 -
1410681.345880*(2-5)^2 -

觀察後面的(i-j)^2項。
是不是相當於
\(-0.0625,-0.111111x,-0.25x^2,-1x^3,0x^4,1x^5,0.25x^6,0.111111x^7,0.0625x^8\)
\(*\)
\(4006373.885184,15375036.435759x,1717456.469144x^2,8514941.004912x^3,1410681.345880x^4\)
轉換為卷積就可以了。

Code

#include<cstdio>
#include<cmath>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<complex>
#define debug cout<<"JEFF是我們的紅太陽!!"<<endl;
#define ll long long
using namespace std;
typedef complex<double> cp;
const int N=1e6+5;
const double pi=acos(-1.0);
cp a[N],b[N];
ll l,n,cnt,limit=1,r[N];
int read(){
    int x=0,w=1;char ch=getchar();
    while(ch>'9'||ch<'0'){if(ch=='-')w=-1;ch=getchar();}
    while(ch>='0'&&ch<='9')x=x*10+ch-'0',ch=getchar();
    return x*w;
}

void FFT(cp *A,int type){
    for(int i=0;i<limit;i++)if(i<r[i])swap(A[r[i]],A[i]);
    for(int i=1;i<limit;i<<=1){
        cp wn(cos(pi/i),sin(type*pi/i));
        for(int j=0;j<limit;j+=(i<<1)){
            cp w(1,0),x,y;
            for(int k=0;k<i;k++){
                x=A[k+j];y=A[k+j+i]*w;
                A[k+j]=x+y;A[k+j+i]=x-y;
                w=w*wn;
            }
        }
    }
}

int main(){
    n=read();
    for(int i=0;i<n;i++){
        double x;
        scanf("%lf",&x);
        b[i]=x;
    }
    for(int i=n-1;i>=1;i--){
        a[cnt]=-(1.0/(1.0*i*i));cnt++;
    }a[cnt]=0;cnt++;
    for(int i=1;i<=n-1;i++){
        a[cnt]=(1.0/(1.0*i*i));cnt++;
    }cnt--;
    while(cnt+n>=limit)limit<<=1,l++;
    for(int i=1;i<limit;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    FFT(a,1);FFT(b,1);
    for(int i=0;i<limit;i++)a[i]*=b[i];
    FFT(a,-1);for(int i=0;i<limit;i++)a[i]/=limit;
    for(int i=n-1;i<=2*n-2;i++)printf("%.3lf\n",(double)(a[i].real()));
    return 0;
}