FFT的应用
2024-09-06 07:21:14
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\}
\]其中\(1\leq n \leq 5\times 10^4,1\leq a_i,b_i\leq 100\)
题解：
首先，有一个结论：两个手环增加非负整数亮度，等于其中一个增加一个整数亮度（可以为负）
我们令增加量为\(x\)，旋转以后的原数列为\(\{a\}\{b\}\)那么现在的费用就是：
\[\sum_{i=1}^n\left(a_i+x-b_i\right)^2
\]
\]
我们把第i项拿出来拆开，得到：
\[\left(a_i+x-b_i\right)^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+2x\left(\sum_{i=1}^na_i-\sum_{i=1}^nb_i\right)-2\sum_{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\times2\)项的最大值，就是
\[\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\)。
\(1\leq n\leq 2 \times 10^5, -10^9\leq 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_{0\leq i < |s|}t_{i+k}(s_i-t_{i+k})^2=0
\]
\]
结合计算卷积，很容易算出左式，即可完成匹配
模板：P4173 残缺的字符串
要求匹配两个字符串A,B，两者都有通配符
\[\sum_{0\leq 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\n",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\)。\(1\leq |s|,|t|, k\leq 2\times 10^5\)，字符集大小为\(4\)。
题解
考虑分别处理每种字符。处理字符\(c\) 时，将\(t\) 中所有和字符为\(c\)的位置距离不超过\(d\) 的位置设为\(1\)，其余位置设为\(0\)；将\(s\) 中所有字符为\(c\) 的位置设为\(1\)，其余位置设为\(0\)。那么在\(k\) 处字符\(c\)的匹配位置个数为
\[\sum_{0\leq i\leq|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;
}
最新文章
- osg渲染数据高程文件
- [zz] be similar with和be similar to的区别
- HTML head 头标签
- TestNG学习-002-annotaton 注解概述及其执行顺序
- docker无法连接进程
- deque双向队列（转）
- float编码杂谈
- 王垠：Lisp 已死，Lisp 万岁！
- AngularJS - 依赖注入(Dependency Injection)
- zz:linux下rz,sz的安装方法
- SpringCloud入门之YAML格式文件规范学习
- IDE安装Lombok插件提高开发效率
- asp.net 微信JsSDK
- linux中一些简便的命令之cut
- c#词频统计命令行程序
- AR介绍
- eclipse .properties插件
- Alpha Version Release Of Teamwork: Appendix 1 BUG BASH
- [转]Linux常用命令大全
- BackBone.js之Router