[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
",(double)(a[i].real()));
    return 0;
}