FFT的应用
——讲稿
概述
FFT的模板很简单，大家都会背，于是出题的空间就在于建模了。FFT的题目难在建模，往往需要将问题抽象出来，经过一系列转化后得到乘积式的和，再赋予式子各个项的系数一定的意义即可。
一些链接：
基本形式
对于类似(displaystyle sum_{i+j=N+k}a_ib_j)的式子，可以直接通过FFT计算。
其中N是定值，表示元素个数；k是变量，是题目中的系数，a和b是间接已知数组，当a与b的系数和为定值时，可将一个翻转，否则直接计算。
例题：P3723 [AH2017/HNOI2017]礼物
给出循环节长度为(n) 的数组(a, b)，求
[min_{c,k}{sum_{i=1}^n(c+a_{i+k}-b_i)^2} ]其中(1leq n leq 5 imes 10^4,1leq a_i,b_ileq 100)
题解：
首先，有一个结论：两个手环增加非负整数亮度，等于其中一个增加一个整数亮度（可以为负）
我们令增加量为(x)，旋转以后的原数列为({a}{b})那么现在的费用就是：
[sum_{i=1}^nleft(a_i+x-b_i
ight)^2
]
我们把第i项拿出来拆开，得到：
[left(a_i+x-b_i
ight)^2=a_i^2+b_i^2+x^2+2a_ix-2a_ib_i-2b_ix
]
那么原式变成了
[sum_{i=1}^na_i^2+sum_{i=1}^nb_i^2+nx^2+2xleft(sum_{i=1}^na_i-sum_{i=1}^nb_i
ight)-2sum_{i=1}^na_ib_i
]
我们发现，这个式子除了最后一项之外都是确定的,
那么我们只要令最后一项最大，那么就可以得到最小的费用值了现在问题转化为求$$sum_{i=1}^na_ib_i$$
这个形式非常卷积，我们把数列({a})反过来，变成
[sum_{i=1}^na_{n-i+1}b_i
]
这不是一个卷积吗~所以把反过来的数列({a})倍长（为了循环卷积），和数列({b})卷积，得到的项里面的第(n+1)到(n imes2)项的最大值，就是
[sum_{i=1}^na_ib_i
]
的最大值然后把前面的不变项加上，就是答案了
代码
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef long double LD;
const int INF=1e9+7,MAXN=1e6+10;
const LD Pi=3.141592653589793238462643383279;
struct com{
LD real,imag;
com(){
real=imag=0;
}
com(LD x,LD y){
real=x;
imag=y;
}
friend inline com operator+(com x,com y){
return com(x.real+y.real,x.imag+y.imag);
}
friend inline com operator-(com x,com y){
return com(x.real-y.real,x.imag-y.imag);
}
friend inline com operator*(com x,com y){
return com(x.real*y.real-x.imag*y.imag,x.real*y.imag+x.imag*y.real);
}
}f[MAXN],g[MAXN];
int N,M,lim=1,L,rev[MAXN];
inline void FFT(com *a,int type){
for(int i=0;i<lim;i++)
if(i<rev[i])
swap(a[i],a[rev[i]]);
for(int hf=1;hf<lim;hf<<=1){
int len=hf<<1;
com Wn=com(cos(2.0*Pi/len),type*sin(-2.0*Pi/len));
for(int j=0;j<lim;j+=len){
com w=com(1,0);
for(int k=0;k<hf;k++){
com t1=a[j+k],t2=a[j+k+hf]*w;
a[j+k]=t1+t2;
a[j+k+hf]=t1-t2;
w=w*Wn;
}
}
}
if(type==-1)
for(int i=0;i<lim;i++)
a[i].real=a[i].real/lim+0.5;
}
LL a[MAXN],b[MAXN],sum1,ans,sum2;
inline LL func1(LL x){
return N*x*x+2*sum1*x;
}
int main(){
scanf("%d%d",&N,&M);
for(int i=1;i<=N;i++)
scanf("%lld",a+i);
for(int i=1;i<=N;i++)
scanf("%lld",b+i);
for(int i=1;i<=N;i++){
ans+=a[i]*a[i]+b[i]*b[i];
sum1+=a[i]-b[i];
}
ans+=min(func1(floor(-(LD)sum1/N)),func1(ceil(-(LD)sum1/N)));
/*--------------------convolution---------------------*/
for(int i=1;i<=N;i++){
f[i].real=a[N-i+1];
g[i].real=g[i+N].real=b[i];
}
while(lim<=N*3){
lim<<=1;
L++;
}
for(int i=0;i<lim;i++)
rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
FFT(f,1);
FFT(g,1);
for(int i=0;i<lim;i++)
f[i]=f[i]*g[i];
FFT(f,-1);
for(int i=N;i<=2*N;i++)
sum2=max(sum2,(LL)f[i].real);
printf("%lld",ans-2*sum2);
return 0;
}
直接计算卷积
一些问题经过简单的推导即可推出可以用计算卷积来解决，这类问题多形如对所有不大于(n)的(k)求出某个东西，其中(k)的答案为求某个卷积后结果数组的第(k)项。
例题：CF993E
给出一个大小为(n)的数组(a)和一个数(x)，对于(0)和(n)之间的所有(k)，求有多少个(a)的区间中恰有(k)个数小于(x)。
(1leq nleq 2 imes 10^5, -10^9leq x,ai leq 10^9)
题解
由于(x) 固定，直接把(a)中小于(x) 数的设为(1)，其余数设为(0)，只需要求有多少个区间的和为(k)。设(s) 为(a) 的前缀和，则只需要求有多少个数对(s_i, s_j) 满足(s_j-s_i = k)，或(s_j - k = s)i。设(f_{s_i}) 为(s_i) 的出现次数，则(k) 的答案为
[sum_{i=k}^nf_if_{i-k}=sum_{i=k}^nf_if_{n-i+k}=sum_{i+j=n+k}f_if_{n-j}
]
设(g_i=f_{n-i})，求(f)和(g)的卷积即可。
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int INF=1e9+7,MAXN=1e6+10/*Min:2^20+10*/;
const long double Pi=3.141592653589793238462643383279;
struct COM{
long double real,imag;
COM(){
real=imag=0;
}
COM(long double x,long double y){
real=x;
imag=y;
}
friend COM operator+(COM x,COM y){
return COM(x.real+y.real,x.imag+y.imag);
}
friend COM operator-(COM x,COM y){
return COM(x.real-y.real,x.imag-y.imag);
}
friend COM operator*(COM x,COM y){
return COM(x.real*y.real-x.imag*y.imag,x.real*y.imag+x.imag*y.real);
}
}f[MAXN],g[MAXN];
int N,X,limit=1,lg,rev[MAXN];
inline void FFT(COM *a,int type){
for(int i=0;i<limit;i++)
if(i<rev[i])
swap(a[i],a[rev[i]]);
for(int half=1;half<limit;half<<=1){
int len=half<<1;
COM Wn=COM(cos(Pi/half),sin(Pi/half)*type);
for(int j=0;j<limit;j+=len){
COM w=COM(1,0),x,y;
for(int k=0;k<half;k++){
x=a[j+k],y=a[j+k+half]*w;
a[j+k]=x+y;
a[j+k+half]=x-y;
w=w*Wn;
}
}
}
if(type==-1)
for(int i=0;i<limit;i++)
a[i].real/=(double)limit;
}
LL s[MAXN],ans[MAXN];
char ch;
int main(){
scanf("%d%d",&N,&X);
for(int i=1;i<=N;i++){
scanf("%lld",s+i);
s[i]=(s[i]<X)+s[i-1];
}
for(int i=0;i<N;i++){
f[s[i]].real++;
g[s[N]-s[i+1]].real++;
}
while(limit<=(s[N]<<1)){
limit<<=1;
lg++;
}
for(int i=0;i<limit;i++)
rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg-1));
FFT(f,1);
FFT(g,1);
for(int i=0;i<limit;i++)
f[i]=f[i]*g[i];
FFT(f,-1);
LL cnt=0;
for(int i=1;i<=N;i++)
if(s[i]==s[i-1])
cnt++;
else{
ans[0]+=cnt*(cnt-1)/2+cnt;
cnt=0;
}
ans[0]+=cnt*(cnt-1)/2+cnt;
for(int i=0;i<s[N];i++)
ans[s[N]-i]=f[i].real+1e-2;
for(int i=0;i<=N;i++)
printf("%lld ",ans[i]);
return 0;
}
多项式运算
多项式乘法可以用来表示卷积，而借助多项式的性质，可以分析并解决类型更为广泛的问题。其中，最典型的例子是利用生成函数解决组合计数问题，这往往可以简化推导过程，有时还可以借助专用算法优化复杂度。
字符串匹配
代通配符的字符串匹配
设(s,t)为字符串，其中(t)中某些字符是通配符，可以匹配任意字符，求(s)在(t)中的所有匹配的位置。将通配符设为(0)，其余字符设为非(0)的数，则(s)在(k)处匹配当且仅当
[sum_{0leq i < |s|}t_{i+k}(s_i-t_{i+k})^2=0
]
结合计算卷积，很容易算出左式，即可完成匹配
模板：P4173 残缺的字符串
要求匹配两个字符串A,B，两者都有通配符
[sum_{0leq i < |s|}t_{i+k}s_i(s_i-t_{i+k})^2=0
]
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int INF=1e9+7,MAXN=2e6+100;
const double Pi=3.14159265358979323846;
struct com{
double real,imag;
com(){
real=imag=0;
}
com(double x,double y){
real=x;
imag=y;
}
friend inline com operator+(com x,com y){
return com(x.real+y.real,x.imag+y.imag);
}
friend inline com operator-(com x,com y){
return com(x.real-y.real,x.imag-y.imag);
}
friend inline com operator*(com x,com y){
return com(x.real*y.real-x.imag*y.imag,x.real*y.imag+x.imag*y.real);
}
friend inline com operator*(com x,int y){
return com(x.real*y,x.imag*y);
}
}f[MAXN],g[MAXN],h[MAXN];
int lim=1,L,rev[MAXN];
inline void FFT(com *a,int type){
for(int i=0;i<lim;i++)
if(i<rev[i])
swap(a[i],a[rev[i]]);
for(int hf=1;hf<lim;hf<<=1){
int len=hf<<1;
com Wn=com(cos(Pi/hf),type*sin(Pi/hf));
for(int j=0;j<lim;j+=len){
com w=com(1,0);
for(int k=0;k<hf;k++){
com t1=a[j+k],t2=a[j+k+hf]*w;
a[j+k]=t1+t2;
a[j+k+hf]=t1-t2;
w=w*Wn;
}
}
}
if(type!=1)
for(int i=0;i<lim;i++)
a[i].real=a[i].real/lim;
}
int M,N;
int a[MAXN],b[MAXN],cnt[MAXN];
vector<int> out;
inline void calc(){
while(lim<=N+M){
lim<<=1;
L++;
}
for(int i=0;i<lim;i++)
rev[i]=rev[i>>1]>>1|(i&1)<<(L-1);
for(int i=0;i<M;i++)
f[i].real=a[i]*a[i]*a[i];
for(int i=0;i<N;i++)
g[i].real=b[i];
FFT(f,1); FFT(g,1);
for(int i=0;i<lim;i++)
h[i]=f[i]*g[i];
for(int i=0;i<lim;i++)
f[i]=g[i]=com();
for(int i=0;i<M;i++)
f[i].real=a[i]*a[i];
for(int i=0;i<N;i++)
g[i].real=b[i]*b[i];
FFT(f,1); FFT(g,1);
for(int i=0;i<lim;i++)
h[i]=h[i]-f[i]*g[i]*2;
for(int i=0;i<lim;i++)
f[i]=g[i]=com();
for(int i=0;i<M;i++)
f[i].real=a[i];
for(int i=0;i<N;i++)
g[i].real=b[i]*b[i]*b[i];
FFT(f,1); FFT(g,1);
for(int i=0;i<lim;i++)
h[i]=h[i]+f[i]*g[i];
FFT(h,-1);
for(int i=M-1;i<N;i++)
if(fabs(h[i].real)<0.5)
out.push_back(i-M+2);
}
char str1[MAXN],str2[MAXN];
int main(){
scanf("%d%d %s %s",&M,&N,str1,str2);
for(int i=0;i<M;i++)
a[M-i-1]=str1[i]=='*'?0:str1[i]-'a'+1;
for(int i=0;i<N;i++)
b[i]=str2[i]=='*'?0:str2[i]-'a'+1;
calc();
int siz=out.size();
printf("%d
",siz);
for(int i=0;i<siz;i++)
printf("%d ",out[i]);
return 0;
}
例题：CF528D Fuzzy Search
给出字符串(s,t)和非负整数(d)，求有多少个(k)，满足对于所有(s)的下标(i)，都存在距离(k+i) 不大于(d)的(j)，使得(s_i = t_j)。(1leq |s|,|t|, kleq 2 imes 10^5)，字符集大小为(4)。
题解
考虑分别处理每种字符。处理字符(c) 时，将(t) 中所有和字符为(c)的位置距离不超过(d) 的位置设为(1)，其余位置设为(0)；将(s) 中所有字符为(c) 的位置设为(1)，其余位置设为(0)。那么在(k) 处字符(c)的匹配位置个数为
[sum_{0leq ileq|s|}s_it_{i+k}=sum_{i+j=|s|+k}{t_is_j}
]
可以通过卷积计算。之后只要判断在每个(k) 处的每种字符的匹配位置个数之和是否等于(|s|) 即可。
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef long double LD;
const int INF=1e9+7,MAXN=1e6+10;
const LD Pi=3.141592653589793238462643383279;
struct com{
LD real,imag;
com(){
real=imag=0;
}
com(LD x,LD y){
real=x;
imag=y;
}
friend inline com operator+(com x,com y){
return com(x.real+y.real,x.imag+y.imag);
}
friend inline com operator-(com x,com y){
return com(x.real-y.real,x.imag-y.imag);
}
friend inline com operator*(com x,com y){
return com(x.real*y.real-x.imag*y.imag,x.real*y.imag+x.imag*y.real);
}
}f[MAXN],g[MAXN];
int lim=1,L,rev[MAXN];
inline void FFT(com *a,int type){
for(int i=0;i<lim;i++)
if(i<rev[i])
swap(a[i],a[rev[i]]);
for(int hf=1;hf<lim;hf<<=1){
int len=hf<<1;
com Wn=com(cos(2.0*Pi/len),type*sin(-2.0*Pi/len));
for(int j=0;j<lim;j+=len){
com w=com(1,0);
for(int k=0;k<hf;k++){
com t1=a[j+k],t2=a[j+k+hf]*w;
a[j+k]=t1+t2;
a[j+k+hf]=t1-t2;
w=w*Wn;
}
}
}
if(type==-1)
for(int i=0;i<lim;i++)
a[i].real=a[i].real/lim+0.5;
}
int N,M,K,ans,cnt[MAXN];
char s[MAXN],t[MAXN];
inline void calc(char ch){
memset(f,0,sizeof(f));
memset(g,0,sizeof(g));
for(int i=0,l=-INF;i<N;i++){
if(s[i]==ch)
l=i;
if(i-l<=K)
f[i].real=1;
}
for(int i=N-1,l=INF;i>=0;i--){
if(s[i]==ch)
l=i;
if(l-i<=K)
f[i].real=1;
}
for(int i=0;i<M;i++)
g[i].real=t[M-i-1]==ch;
FFT(f,1);
FFT(g,1);
for(int i=0;i<lim;i++)
f[i]=f[i]*g[i];
FFT(f,-1);
for(int i=0;i<N;i++)
cnt[i]+=floor(f[i+M-1].real);
}
int main(){
scanf("%d%d%d %s %s",&N,&M,&K,s,t);
while(lim<=N+M-2){
lim<<=1;
L++;
}
for(int i=0;i<lim;i++)
rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
for(int i=0;i<4;i++)
calc("AGCT"[i]);
for(int i=0;i<N;i++)
ans+=cnt[i]>=M;
printf("%d",ans);
return 0;
}