FFT入门及其应用

系列博客 https://www.cnblogs.com/2016gdgzoi471/p/9476883.html 

FFT原理详解，包含大数乘法的代码

https://www.cnblogs.com/RabbitHu/p/FFT.html

多项式乘法：就是套一套板子：

#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
#define ll long long
#define N 2000005

struct complex{
    double x,y;
    complex (double xx=0,double yy=0){x=xx,y=yy;}
};
complex operator + (complex a,complex b){ return complex(a.x+b.x , a.y+b.y);}
complex operator - (complex a,complex b){ return complex(a.x-b.x , a.y-b.y);}
complex operator * (complex a,complex b){ return complex(a.x*b.x-a.y*b.y , a.x*b.y+a.y*b.x);}

typedef complex cp;
const double pi = acos(-1.0);
inline int read(){
    char c=getchar();int x=0,f=1;
    while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
    return x*f;
}
ll A[N],B[N],lena,lenb,n,res[N];
cp a[N],b[N],omg[N],inv[N];

void init(){//预处理出单位根及其倒数 
    for(register int i=0;i<n;i++)
        omg[i]=cp(cos(2*pi/n*i),sin(2*pi/n*i));
    for(register int i=0;i<n;i++)
        inv[i]=cp(cos(2*pi/n*i),-sin(2*pi/n*i));
}

void fft(cp *a, cp *omg){
    int lim = 0;
    while((1<<lim)<n)lim++;    //向下分治的次数 
    for(register int i=0;i<n;i++){    //先把系数移到最后一层递归的位置 
        int t=0;
        for(register int j=0;j<lim;j++)
            if((i>>j) & 1) t|=(1<<(lim-j-1));
        if(i<t)swap(a[i],a[t]);
    }
    for(register int l=2;l<=n;l*=2){    //每个递归区间的长度 
        int m=l/2;
        for(cp *p=a;p!=a+n;p+=l)
            for(register int i=0;i<m;i++){//蝴蝶操作 
                cp t=omg[n/l*i]*p[i+m];
                p[i+m]=p[i]-t;
                p[i]=p[i]+t; 
            }
    }
    
}

//a0,a1,a2,a3...an

int main(){
    lena=read();lenb=read();
    lena++;lenb++;
    n=1;while(n<lena+lenb)n<<=1;
    for(register int i=0;i<lena;i++)A[i]=read();
    for(register int i=0;i<lenb;i++)B[i]=read();
    for(register int i=0;i<lena;i++)a[i]=cp(A[i],0);
    for(register int i=0;i<lenb;i++)b[i]=cp(B[i],0);
    init();
    fft(a,omg);fft(b,omg);
    for(register int i=0;i<n;i++)a[i]=a[i]*b[i];
    fft(a,inv);
        
    for(register int i=0;i<n;i++)
        res[i]=floor(a[i].x/n+0.5);
    for(register int i=0;i<lena+lenb-1;i++)
        cout<<res[i]<<" ";
        
}

多项式乘法的卷积形式

https://www.cnblogs.com/RabbitHu/p/BZOJ2194.html

FFT在字符串相似度匹配上的应用：

https://blog.csdn.net/qq_31759205/article/details/80060918

例题：

GYM101667H 题解： https://blog.csdn.net/ACTerminate/article/details/79325574

#include<bits/stdc++.h>
#include<complex>
using namespace std;
#define ll long long 
#define N 400005
typedef complex<double> cp;

int lena,lenb;
char A[N],B[N];

cp a[N],b[N],omg[N],inv[N];
int n,sum[N],F[N];
const double pi = acos(-1.0);
void init(){
    for(int i=0;i<n;i++){
        omg[i]=cp(cos(2*pi/n*i),sin(2*pi/n*i));
        inv[i]=conj(omg[i]);
    }
}
void fft(cp *a, cp * omg){
    int lim=0;
    while((1<<lim)<n)lim++;
    for(int i=0;i<n;i++){
        int t=0;
        for(int j=0;j<lim;j++)
            if((i>>j) & 1) t|=(1<<(lim-j-1));
        if(i<t)swap(a[i],a[t]);    
    }
    for(int l=2;l<=n;l*=2){
        int m=l/2;
        for(cp *p=a;p!=a+n;p+=l){
            for(int i=0;i<m;i++){
                cp t=omg[n/l*i]*p[i+m];
                p[i+m]=p[i]-t;
                p[i]=p[i]+t;
            }
        }
    }
}
    
void calc(char ch){
    memset(a,0,sizeof a);
    memset(b,0,sizeof b);
    for(int i=0;i<lena;i++)
        if(A[i]==ch) a[i].real(1);
        else a[i].real(0);
    for(int i=0;i<lenb;i++)
        if(B[i]==ch) b[i].real(1);
        else b[i].real(0);
    fft(a,omg);fft(b,omg);
    for(int i=0;i<n;i++)
        a[i]=a[i]*b[i];
    fft(a,inv);
    for(int i=0;i<n;i++)
        sum[i]=floor(a[i].real()/n+0.5);
    for(int i=0;i<lena+lenb-1;i++)
        F[i]+=sum[i];
} 

int main(){
    cin>>lena>>lenb;
    n=1;
    while(n<lena+lenb)n<<=1;
    init();
    scanf("%s%s",A,B);
    
    for(int i=0;i<lena;i++){
        if(A[i]=='S')A[i]='R';
        else if(A[i]=='R')A[i]='P';
        else if(A[i]=='P')A[i]='S';
    }
    for(int i=0,j=lenb-1;i<j;i++,j--)
        swap(B[i],B[j]);
        
    //对每种字符单独求，再加起来 
    calc('R');calc('P');calc('S');
    
    //F[k]=sum{A[i]*B[k-i]},表示从A[k-(lenb-1)]开始匹配的相似度 
    //因为是匹配一整个B串，所以k>=lenb-1 
    int Max=0;
    for(int i=lenb-1;i<lena+lenb-1;i++)
        Max=max(Max,F[i]);
    cout<<Max<<'
';
} 

2018宁夏网络赛H题：https://nanti.jisuanke.com/t/A1676

题解：https://blog.csdn.net/qq_38891827/article/details/80281151

一份不错的板子：

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#include<complex>
using namespace std;
#define maxn 2000005

const double pi=acos(-1.0);

typedef complex<double>cp;
cp a[maxn],b[maxn];

int S,T,n,m,L,R[maxn],ans[maxn];
long long F[maxn];
char s[maxn],t[maxn];
double f[maxn],g[maxn];

void FFT(cp a[maxn],int opt)
{
    int lim = 0;
    while((1 << lim) < n) lim++;
    for(register int i = 0; i < n; i++){
        int t = 0;
        for(register int j = 0; j < lim; j++)
            if((i >> j) & 1) t |= (1 << (lim - j - 1));
        if(i < t) swap(a[i], a[t]); // i < t 的限制使得每对点只被交换一次（否则交换两次相当于没交换）
    }
    for (int k=1;k<n;k<<=1)
    {
        cp wn=cp(cos(pi/k),opt*sin(pi/k));
        for (int i=0;i<n;i+=(k<<1))
        {
            cp w=cp(1,0);
            for (int j=0;j<k;++j,w=w*wn)
            {
                cp x=a[i+j],y=w*a[i+j+k];
                a[i+j]=x+y,a[i+j+k]=x-y;
            }
        }
    }
}
void calc(char ch1,char ch2,char ch3)
{
    memset(a,0,sizeof a);
    memset(b,0,sizeof b);
    for(int i=0;i<T;i++)
        if(t[i]==ch2 || t[i]==ch3)a[i].real(1);
    for(int i=0;i<S;i++)
        if(s[i]==ch1)b[i].real(1);
    
    FFT(a,1);FFT(b,1);
    for (int i=0;i<=n;++i) a[i]=a[i]*b[i];
    FFT(a,-1);
    for (int i=S-1;i<T;++i)
    {
         F[i]+=(long long)(a[i].real()/n+0.5);//对于每种匹配位置累计赢的次数
    }
}
int main()
{
    scanf("%s%s",t,s);//S短,T长
    S=strlen(s);
    T=strlen(t);
    for (int i=0;i<S/2;++i) swap(s[i],s[S-i-1]);//短串逆置

    m=S+T-2;
    for (n=1;n<=m;n<<=1) ++L;

        
    calc('S','P','L');
    calc('P','R','K');
    calc('R','L','S');
    calc('L','K','P');
    calc('K','S','R');

    long long  ans=0;
    for (int i=S-1;i<T;++i)//这个范围自己考虑一下就好了
      ans=max(ans,F[i]);//所有位置取max
    printf("%lld
",ans);
}

 bzoj4259:同样是逆置字符串的思路，还多了一个分段求FFT的思路

但是bzoj上不能交了https://www.cnblogs.com/clrs97/p/4814499.html

16香港网络赛：FFT求A[i]+A[j]=A[k]的对数，FFT+容斥+坐标偏移

题解：https://blog.csdn.net/qq_38891827/article/details/80281151

/*
F[k]=sum{a[i]*b[k-i]}
由于本题要把坐标向右偏移LOW，得
    F[k+2*LOW]=sum{a[i+LOW]*b[k-i+LOW]}

枚举每个数A[i]    
    那么F[A[i]+2*LOW]就表示x^A[i]出现的次数 
再容斥掉一些答案：
    1.自身*自身：F[(A[i]+LOW)*2]--;
    2.A[i]!=0，每当A[j]=0时：
        所有(i,j,i),(j,i,i)都会被算贡献，所以减掉0的个数 
    3.A[i]=0，每当A[j]=0的情况：
        所有(i,j,i),(j,i,i)都会被算贡献，所以减掉0的个数-1 
*/
#include<bits/stdc++.h>
#include<complex>
using namespace std;
#define ll long long 
#define N 1000005
#define LOW 50000

ll F[N],m,A[N],cnt[N],n;
typedef complex<double> cp;
cp a[N],b[N];
const double pi = acos(-1.0);

void FFT(cp a[N], int opt){
    int lim=0;
    while((1<<lim)<n)lim++;
    for(int i=0;i<n;i++){
        int t=0;
        for(int j=0;j<lim;j++)
            if((i>>j) & 1)t|=(1<<(lim-j-1));
        if(i<t)swap(a[i],a[t]);
    }
    for(int k=1;k<n;k<<=1){
        cp wn=cp(cos(pi/k),opt*sin(pi/k));//原根 
        for(int i=0;i<n;i+=(k<<1)){
            cp w=cp(1,0);
            for(int j=0;j<k;j++,w=w*wn){
                cp x=a[i+j],y=w*a[i+j+k];
                a[i+j]=x+y;
                a[i+j+k]=x-y;
            }
        }
    }
}

int zero;
int main(){
    cin>>m;
    for(int i=1;i<=m;i++){
        scanf("%lld",&A[i]);
        cnt[A[i]+LOW]++;
        if(A[i]==0)zero++;
    }
    n=1;
    while(n<200000)n<<=1;
    for(int i=0;i<n;i++)
        a[i].real(cnt[i]),b[i].real(cnt[i]);
        
    FFT(a,1);FFT(b,1);
    for(int i=0;i<n;i++)a[i]*=b[i];
    FFT(a,-1);
    
    for(int i=0;i<n;i++)F[i]=(ll)(a[i].real()/n+0.5);
    for(int i=1;i<=m;i++)F[(A[i]+LOW)*2]--;
    ll ans=0;
    for(int i=1;i<=m;i++){//每个数出现的次数 
        ans+=F[A[i]+LOW*2];
        //减去0的影响 
        if(A[i]==0)
            ans-=(zero-1)*2;
        else ans-=zero*2;
    }
    
    cout<<ans<<'
';
}