【Bzoj3527】力

来源和评测点

[Zjoi2014]
Bzoj3527
Caioj1451
Caioj1451

题目

【题目大意】
给出n个数qi,给出Fj的定义如下:

令Ei=Fi/qi,求Ei
【输入格式】
第一行一个整数n
接下来n行每行输入一个数,第i行表示qi
n≤100000
0<qi<1000000000
【输出格式】
n行,第i行输出Ei。与标准答案误差不超过1e-2即可。
【输入样例1】
3
124.44513
547.26513
5789.147852
【输出样例1】
-1994.552
-5664.703
578.376
【输入样例2】
5
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880
【输出样例2】
-16838672.693
3439.793
7509018.566
4595686.886
10903040.872

分析

快速傅里叶变换教程和题目分类参见:
【OI之路】11更高级数论-3快速傅里叶变换
其他快速傅里叶变换题目参见:Tag-快速傅里叶变换

这道题还是挺有难度的
基本思路:逐步化简公式转化为卷积的形式从而能用FFT解决
下标统一从0开始,翻转时也会用到

$$
题目要求:
F_i=
\sum_{i < k} \frac{q_i q_k}{(i-k)^2}-
\sum_{i > k} \frac{q_i q_k}{(k-i)^2}
$$

$$
而E_i=\frac{F_i}{q_i}
$$

$$
显而易见
E_i=
\sum_{i < k} \frac{q_k}{(i-k)^2}-
\sum_{i > k} \frac{q_k}{(k-i)^2}
$$

$$
明确界限
E_i=
\sum_{k=0}^{i-1} \frac{q_k}{(i-k)^2}-
\sum_{k=i+1}^{n-1} \frac{q_k}{(k-i)^2}
$$

$$
设A_t=q_t,B_t=\frac{1}{t^2}(精髓1)
$$

$$
E_i=
\sum_{k=0}^{i-1} A_k B_{i-k}-
\sum_{k=i+1}^{n-1} A_k B_{k-i}
$$

$$
设B_0=0(小细节),这样前面的式子就和卷积一样了,重点放后面
$$

$$
E_i=
\sum_{k=0}^{i} A_k B_{i-k}-
\sum_{k=i}^{n-1} A_k B_{k-i}
$$

考虑$A’t=A{t+i}$,$B’_{(n-i-1)-k}=B_k$(精髓2,一个前移,一个翻转)

$$
E_i=
\sum_{k=0}^{i} A_k B_{i-k}-
\sum_{k=0}^{n-i-1} A’k B’{(n-i-1)-i}
$$

$$
E_i=FFT1(i)-FFT2(n-i-1)
$$
其中两个FFT只是A不同

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
//Zory-2017
//*******************头文件*******************
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
//*******************全局常量*******************
const double PI=acos(-1.0);
const int MAXN=4*100000;
//*******************全局定义*******************
struct complex
{
double r,i;//real,imag
complex() {r=i=0;}
complex(double rr,double ii) {r=rr;i=ii;}
complex operator + (complex b) {return complex(r+b.r,i+b.i);}
complex operator - (complex b) {return complex(r-b.r,i-b.i);}
complex operator * (complex b) {return complex(r*b.r-i*b.i,r*b.i+i*b.r);}
//(a+bi)(c+di)=ac+adi+bci+bdi^2=(ac-bd)+(ad+bc)i
//baike.baidu.com/item/复数运算法则
};
//*******************实现*******************
complex omega[MAXN+10];
void CalcOmega(int n,int op)
{
double t=2*PI*op/n;
omega[0]=complex(1,0);
omega[1]=complex(cos(t),sin(t));
for(int k=2;k<=n/2-1;k++) omega[k]=omega[k-1]*omega[1];
}
int R[MAXN+10];
void fft(complex *a,int n,int op)
{
for(int i=0;i<=n-1;i++) if(i<R[i]) swap(a[i],a[R[i]]);
//注意判断大小不然又换回来了
for(int L=1;L<=n/2;L*=2)//合并前长度
{
CalcOmega(L*2,op);
for(int j=0;j<=n-1;j+=L*2)//开头,长度=L*2
{
for(int k=0;k<=L-1;k++)
{
complex x=a[j+k];
complex y=omega[k]*a[j+L+k];
a[j+k]=x+y;
a[j+L+k]=x-y;
//蝴蝶操作
}
}
}
if(op==-1) for(int i=0;i<=n-1;i++) a[i].r=a[i].r/double(n);
}
//*******************接口*******************
void FFT(complex *a,complex *b,complex *c,int n)
{
fft(a,n,1);fft(b,n,1);//求值
for(int i=0;i<=n-1;i++) c[i]=a[i]*b[i];//点值乘法
fft(c,n,-1);//插值
}
//*******************主函数*******************
complex a[MAXN+10],b[MAXN+10];
complex f1[MAXN+10],f2[MAXN+10];
double q[MAXN+10];
int main()
{
int n;scanf("%d",&n);
for(int i=0;i<=n-1;i++) scanf("%lf",&q[i]);
int fn=n-1,fm=n-1;
int sn=1,lg=0;
while(sn<fn+fm+1) sn*=2,lg++;
for(int i=0;i<=sn-1;i++) R[i]=( R[i>>1]>>1 ) | ( (i&1)<<(lg-1) );
//nlogn翻转二进制
for(int i=0;i<=n-1;i++) a[i].r=q[i];b[0].r=0;
for(int i=1;i<=n-1;i++) b[i].r=1.0/double(i)/double(i);
FFT(a,b,f1,sn);
for(int i=0;i<=sn-1;i++) a[i]=complex(),b[i]=complex();
for(int i=0;i<=n-1;i++) a[i].r=q[n-i-1];b[0].r=0;
for(int i=1;i<=n-1;i++) b[i].r=1.0/double(i)/double(i);
FFT(a,b,f2,sn);
for(int i=0;i<=n-1;i++) printf("%.3lf\n",f1[i].r-f2[n-i-1].r);
}

本文基于 知识共享署名-相同方式共享 4.0 国际许可协议发布
本文地址:http://zory.cf/2017-12/力.html
转载请注明出处,谢谢!

哪怕是一杯奶茶,也将鼓励我继续创作!
0%