FFT的应用

FFT的应用

——讲稿

概述

FFT的模板很简单，大家都会背，于是出题的空间就在于建模了。FFT的题目难在建模，往往需要将问题抽象出来，经过一系列转化后得到乘积式的和，再赋予式子各个项的系数一定的意义即可。

一些链接：

FFT快速傅里叶变换

NTT数论变换

MTT：任意模数NTT

分治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;
}

给出字符串(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;
}
【推广】 免费学中医，健康全家人
原文地址：https://www.cnblogs.com/guoshaoyang/p/11296027.html