[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;
}
上一篇:在某些情况下明明添加了引用,为何VS还报错"XXX"不存在类型或命名空间(是否缺少程序集引用)


下一篇:如何获取xml循环体报文内容