FFT及相关例题

之前学过FFT,但是感觉自己仅仅停留在简单的多项式直接卷积,于是最近找了点题来做,在这里做一下总结。

FFT简单讲解

入门讲解似乎非常麻烦……我就之讲讲实现上的细节吧……

  • IDFT完之后一定要记得除$n$……因为这个错了好几次……
  • FFT数组最好开到4倍……

例题

这三道题都算是自己做出来的,自我感觉不错?

BZOJ3527 [ZJOI2014]力

题目链接:http://www.lydsy.com/JudgeOnline/problem.php?id=3527

题意

给出$n$个数$q_i$,给出$F_j$的定义如下
$$F_j=\sum_{i<j}\frac{q_iq_j}{(i-j)^2}-\sum_{j<i}\frac{q_iq_j}{(i-j)^2}$$
令$E_i=\frac{F_i}{q_i}$,求$E_i$。
$n \leqslant 1e5$。

题解

首先把那个$q_i$除掉,然后其实就是卷积的形式了(雾)。把$q_i$视作数列{$a_i$},数列{$b_i$}=$\frac{1}{i^2}$,那么可以发现前一个和式就可以用$a_i$与$b_i$卷积。而后一个式子只需把$b_i$倒过来就可以了。
时间复杂度$O(nlogn)$。

代码

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
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long ll;
const int INF=1e9;
const long double eps=1e-9;
const long double Pi=acos(-1.0);
const int maxn=4e5+10;
struct Complex{
long double re,im;
};
Complex operator + (const Complex &lhs,const Complex &rhs){
return (Complex){lhs.re+rhs.re,lhs.im+rhs.im};
}
Complex operator - (const Complex &lhs,const Complex &rhs){
return (Complex){lhs.re-rhs.re,lhs.im-rhs.im};
}
Complex operator * (const Complex &lhs,const Complex &rhs){
return (Complex){lhs.re*rhs.re-lhs.im*rhs.im,lhs.re*rhs.im+lhs.im*rhs.re};
}
struct Complex a[maxn],b[maxn];
long double q[maxn],ans[maxn];
int r[maxn];
int n;
inline int read(){
int x=0,flag=1;
char ch=getchar();
while(!isdigit(ch) && ch!='-')ch=getchar();
if(ch=='-')flag=-1,ch=getchar();
while(isdigit(ch))x=(x<<3)+(x<<1)+ch-'0',ch=getchar();
return x*flag;
}
inline void FFT(Complex P[],int opt){
int i,j,k;
for(i=0;i<n;i++)if(i<r[i])swap(P[i],P[r[i]]);
for(i=2;i<=n;i<<=1){
Complex Wi=(Complex){cos(2*Pi/i),opt*sin(2*Pi/i)};
int p=i/2;
for(j=0;j<n;j+=i){
Complex x=(Complex){1,0};
for(k=0;k<p;k++,x=x*Wi){
Complex u=P[j+k],v=x*P[j+k+p];
P[j+k]=u+v;
P[j+k+p]=u-v;
}
}
}
}
int main(){
int i,j,k,m,tmp;
#ifndef ONLINE_JUDGE
freopen("BZOJ3527.in","r",stdin);
freopen("BZOJ3527.out","w",stdout);
#endif
n=read();tmp=n;m=2*n-1;
for(i=1;i<=n;i++)scanf("%Lf",&q[i]),a[i-1].re=q[i];
for(i=1;i<n;i++)b[i-1].re=1.0/(1ll*(n-i)*(n-i));
int cnt=0;
for(n=1;n<=m;n<<=1)cnt++;
for(i=0;i<n;i++)r[i]=(r[i>>1]>>1) | (i & 1) * (1<<(cnt-1));
FFT(a,1);FFT(b,1);
for(i=0;i<n;i++)a[i]=a[i]*b[i];
FFT(a,-1);
swap(n,tmp);
for(i=0;i<n-1;i++)ans[i]-=a[i+n-1].re/tmp;
memset(a,0,sizeof(a));memset(b,0,sizeof(b));
for(i=1;i<=n;i++)a[i-1].re=q[i];
for(i=1;i<n;i++)b[i-1].re=1.0/(1ll*i*i);
swap(n,tmp);
FFT(a,1);FFT(b,1);
for(i=0;i<n;i++)a[i]=a[i]*b[i];
FFT(a,-1);
for(i=1;i<tmp;i++)ans[i]+=a[i-1].re/n;
for(i=0;i<tmp;i++)printf("%.6Lf\n",ans[i]);
return 0;
}

BZOJ4827 [HNOI2017]礼物

题目链接:http://www.lydsy.com/JudgeOnline/problem.php?id=4827

题意

我的室友最近喜欢上了一个可爱的小女生。马上就要到她的生日了,他决定买一对情侣手环,一个留给自己,一个送给她。每个手环上各有$n$个装饰物,并且每个装饰物都有一定的亮度。但是在她生日的前一天,我的室友突然发现他好像拿错了一个手环,而且已经没时间去更换它了!他只能使用一种特殊的方法,将其中一个手环中所有装饰物的亮度增加一个相同的自然数$c$(即非负整数)。并且由于这个手环是一个圆,可以以任意的角度旋转它,但是由于上面 装饰物的方向是固定的,所以手环不能翻转。需要在经过亮度改造和旋转之后,使得两个手环的差异值最小。在将两个手环旋转且装饰物对齐了之后,从对齐的某个位置开始逆时针方向对装饰物编号$1,2,…,n$,其中$n$为每个手环的装饰物个数,第$1$个手环的$i$号位置装饰物亮度为$x_i$,第$2$个手环的$i$号位置装饰物亮度为$y_i$,两个手环之间的差异值为: $\sum_{i=1}^{n}(x_i-y_i)^2$。麻烦你帮他计算一下,进行调整(亮度改造和旋转),使得两个手环之间的差异值最小, 这个最小值是多少呢?注意在将手环改造之后,装饰物的亮度可以大于$m$。
$1 \leqslant n \leqslant 50000,1 \leqslant m \leqslant 100,a_i \leqslant m$。

题解

首先要解决同时加一个值的问题。不难发现,如果当两个手环的亮度和之差最小时,两个手环的差异值最小。具体证明用二次函数不难得到,在此略去。我们观察差异值的计算式:$\sum_{i=1}^{n}(x_i-y_i)^2$。当$\sum_{i=1}^{n}x_i$与$\sum_{i=1}^{n}y_i$的值已经确定时,为了最小化差异值,我们要最大化$\sum_{i=1}^{n}x_i \times y_i$。这个形式有点像卷积的形式,于是我们考虑FFT。不难发现,当我们把$y$序列反转过来时,这个式子就是卷积的形式了,并且同时恰好在一次FFT中解决了所有的旋转的问题。于是问题得到解决。
时间复杂度$O(nlogn)$。

代码

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
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long ll;
const int INF=2e9;
const long double eps=1e-9;
const double Pi=acos(-1.0);
const int maxn=2e5+10;
struct Complex{
double re,im;
};
Complex operator + (const Complex &lhs,const Complex &rhs){
return (Complex){lhs.re+rhs.re,lhs.im+rhs.im};
}
Complex operator - (const Complex &lhs,const Complex &rhs){
return (Complex){lhs.re-rhs.re,lhs.im-rhs.im};
}
Complex operator * (const Complex &lhs,const Complex &rhs){
return (Complex){lhs.re*rhs.re-lhs.im*rhs.im,lhs.re*rhs.im+lhs.im*rhs.re};
}
struct Complex A[maxn],B[maxn];
int a[maxn],b[maxn],r[maxn],ans[maxn];
int n;
inline int read(){
int x=0,flag=1;
char ch=getchar();
while(!isdigit(ch) && ch!='-')ch=getchar();
if(ch=='-')flag=-1,ch=getchar();
while(isdigit(ch))x=(x<<3)+(x<<1)+ch-'0',ch=getchar();
return x*flag;
}
inline void FFT(Complex P[],int opt){
int i,j,k;
for(i=0;i<n;i++)if(i<r[i])swap(P[i],P[r[i]]);
for(i=2;i<=n;i<<=1){
Complex Wi=(Complex){cos(2*Pi/i),opt*sin(2*Pi/i)};
int p=i/2;
for(j=0;j<n;j+=i){
Complex x=(Complex){1,0};
for(k=0;k<p;k++,x=x*Wi){
Complex u=P[j+k],v=x*P[j+k+p];
P[j+k]=u+v;
P[j+k+p]=u-v;
}
}
}
}
int main(){
int i,j,k,m,tmp;
#ifndef ONLINE_JUDGE
freopen("BZOJ4827.in","r",stdin);
freopen("BZOJ4827.out","w",stdout);
#endif
n=read();m=read();
tmp=n;
int finalans=0,sum1=0,sum2=0,flag=1;
for(i=1;i<=n;i++)a[i]=read(),sum1+=a[i];
for(i=1;i<=n;i++)b[i]=read(),sum2+=b[i];
if(sum1>sum2)swap(sum1,sum2),flag=-1;
int delta=(sum2-sum1)/n+((sum2-sum1)%n>n/2);
delta*=flag;
if(delta>=0)for(i=1;i<=n;i++)a[i]+=delta;
else for(i=1;i<=n;i++)b[i]-=delta;
for(i=1;i<=n;i++)A[i-1].re=a[i],B[i-1].re=b[n-i+1],finalans+=a[i]*a[i]+b[i]*b[i];
int cnt=0;k=n+n;
for(n=1;n<=k;n<<=1)cnt++;
for(i=1;i<=n;i++)r[i]=(r[i>>1]>>1) | ((i & 1) * (1<<(cnt-1)));
FFT(A,1);FFT(B,1);
for(i=0;i<n;i++)A[i]=A[i]*B[i];
FFT(A,-1);
for(i=0;i<n;i++)ans[i]=floor(1.0*A[i].re/n+0.5);
swap(tmp,n);
int maxans=0;
for(i=0;i<=n;i++)
maxans=max(maxans,ans[i]+ans[i+n]);
printf("%d\n",finalans-2*maxans);
return 0;
}

BZOJ3160 万径人踪灭

题目链接:http://www.lydsy.com/JudgeOnline/problem.php?id=3160

题意

给定一个长度为$n$且只含ab的字符串。求满足以下条件的子序列的个数:

  • 这个子序列是一个回文串
  • 这个子序列的对应的位置是对称的
  • 这个子序列不是连续的

方案数对$1e9+7$取模。
$n \leqslant 100000$。

题解

对于第三个条件,我们不难得到转化:跑一遍manacher即可求出以每个字符串为中心的回文串长度,进而可以求出连续的子序列且满足前两个条件的方案数。而对于前两个条件,我们如果求出了求出关于每个位置对称且字符相同的位置的对数,即可轻松计算出最终答案。而对于这个问题,我们分别令a=0,b=1a=1,b=0并各自与自己翻转过来的序列卷积一下,即可分别求出关于每个位置对称且字符为a与为b的位置的对数。问题得到解决。
时间复杂度$O(nlogn)$。

代码

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
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long ll;
const int INF=1e9;
const int mod=1e9+7;
const long double eps=1e-9;
const long double Pi=acos(-1.0);
const int maxn=4e5+10;
struct Complex{
long double re,im;
};
Complex operator + (const Complex &lhs,const Complex &rhs){
return (Complex){lhs.re+rhs.re,lhs.im+rhs.im};
}
Complex operator - (const Complex &lhs,const Complex &rhs){
return (Complex){lhs.re-rhs.re,lhs.im-rhs.im};
}
Complex operator * (const Complex &lhs,const Complex &rhs){
return (Complex){lhs.re*rhs.re-lhs.im*rhs.im,lhs.re*rhs.im+lhs.im*rhs.re};
}
struct Complex a[maxn],b[maxn];
char str[maxn],s[maxn];
int n,P[maxn],ans[maxn],r[maxn];
inline int read(){
int x=0,flag=1;
char ch=getchar();
while(!isdigit(ch) && ch!='-')ch=getchar();
if(ch=='-')flag=-1,ch=getchar();
while(isdigit(ch))x=(x<<3)+(x<<1)+ch-'0',ch=getchar();
return x*flag;
}
inline ll qpow(ll base,ll power){
ll pans=1;
while(power){
if(power & 1ll)pans=pans*base%mod;
base=base*base%mod;
power/=2ll;
}
return pans;
}
inline void FFT(Complex A[],int opt){
int i,j,k;
for(i=0;i<n;i++)if(i<r[i])swap(A[i],A[r[i]]);
for(i=2;i<=n;i<<=1){
Complex Wi=(Complex){cos(2*Pi/i),opt*sin(2*Pi/i)};
int p=i/2;
for(j=0;j<n;j+=i){
Complex x=(Complex){1,0};
for(k=0;k<p;k++,x=x*Wi){
Complex u=A[j+k],v=x*A[j+k+p];
A[j+k]=u+v;
A[j+k+p]=u-v;
}
}
}
}
int main(){
int i,j,k,m,tmp;
#ifndef ONLINE_JUDGE
freopen("BZOJ3160.in","r",stdin);
freopen("BZOJ3160.out","w",stdout);
#endif
scanf("%s",str);
strcpy(s,str);
n=strlen(str);
for(i=n;i>=0;i--){
str[i+i+2]=str[i];
str[i+i+1]='$';
}
str[0]='*';
int id=0;
for(i=2;i<=2*n+1;i++){
if(P[id]+id>i)P[i]=min(P[id*2-i],P[id]+id-i);
else P[i]=1;
while(str[i-P[i]]==str[i+P[i]])P[i]++;
if(P[i]+i>P[id]+id)id=i;
}
for(i=0;i<n;i++)a[i].re=b[i].re=(s[i]=='a'?1:0);
int cnt=0;tmp=n,m=n+n;
for(n=1;n<=m;n<<=1)cnt++;
for(i=0;i<n;i++)r[i]=(r[i>>1]>>1) | (i & 1) * (1<<(cnt-1));
FFT(a,1);FFT(b,1);
for(i=0;i<n;i++)a[i]=a[i]*b[i];
FFT(a,-1);
for(i=0;i<n;i++)ans[i]+=(floor(a[i].re/n+0.5)+1)/2;
memset(a,0,sizeof(a));memset(b,0,sizeof(b));
swap(n,tmp);
for(i=0;i<n;i++)a[i].re=b[i].re=(s[i]=='b'?1:0);
swap(n,tmp);
FFT(a,1);FFT(b,1);
for(i=0;i<n;i++)a[i]=a[i]*b[i];
FFT(a,-1);
for(i=0;i<n;i++)ans[i]+=(floor(a[i].re/n+0.5)+1)/2;
ll finalans=0;
for(i=0;i<n;i++)finalans=(finalans+(qpow(2,ans[i])-1))%mod;
swap(n,tmp);
for(i=0;i<=2*n+1;i++)finalans=(finalans-(P[i]/2))%mod;
printf("%lld\n",(finalans%mod+mod)%mod);
return 0;
}

总结

关于FFT的各类模型,除了可以计算卷积形式的式子之外,还可以处理像例三这样的字符串问题,个人觉得还是很神奇的。BZOJ类似于例三这样的题似乎有很多,但我找到的都是权限题= =。