FFT模板以及简单的FFT入门题总结

66 阅读 0 评论 44 点赞

我是靠谱客的博主温柔金毛，最近开发中收集的这篇文章主要介绍FFT模板以及简单的FFT入门题总结，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

FFT

前言：感谢M哥,Pls,Zjq等等一系列大神的指导，让我大概了解了FFT的基本原理以及能做的事情，下面是FFT的一些简单入门题（我也只会这些题，在这里总结一下。
由于目前题量较少，对于我目前的理解，FFT主要解决两类问题，一类是特殊要求的字符串匹配可构造卷积FFT的，另一类是对两个数组n^2的加法的所有可能结果进行加速。

LOJ 多项式乘法FFT裸题

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
using namespace std;
#define maxn 400005//开四倍所有数组。

const double pi=acos(-1.0);
struct complex
{
    double x,y;
    complex(double X=0,double Y=0)
    {
        x=X,y=Y;
    }
} a[maxn],b[maxn];
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);
}
int S,T,n,m,L,R[maxn],ans[maxn];
long long F[maxn];//存储多项式相乘结果的系数从0-n
double f[maxn],g[maxn];

void FFT(complex a[maxn],int opt)
{
    for (int i=0; i<n; ++i)
        if (i<R[i]) swap(a[i],a[R[i]]);
    for (int k=1; k<n; k<<=1)
    {
        complex wn=complex(cos(pi/k),opt*sin(pi/k));
        for (int i=0; i<n; i+=(k<<1))
        {
            complex w=complex(1,0);
            for (int j=0; j<k; ++j,w=w*wn)
            {
                complex x=a[i+j],y=w*a[i+j+k];
                a[i+j]=x+y,a[i+j+k]=x-y;
            }
        }
    }
}
void calc(int opt)
{
    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=0; i<=n; ++i)//上限为n，在运用时灵活修改。
    {
        F[i]=(long long)(a[i].x/n+0.5)*opt;
        //cout<<F[i]<<endl;
    }
}
int main()
{
    scanf("%d%d",&S,&T);
    for(int i=0; i<=S; i++)
        scanf("%lf",&f[i]);
    for(int i=0; i<=T; i++)
        scanf("%lf",&g[i]);
    m=S+T;
    for (n=1; n<=m; n<<=1) ++L;
    for (int i=0; i<n; ++i)
        R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
    for (int i=0; i<=n; ++i)
        a[i]=complex(1.0*f[i],0.0),b[i]=complex(1.0*g[i],0.0);
    calc(1);
    for(int i=0; i<=S+T; i++)
        printf("%d ",F[i]);
    return 0;
}

FFT第一题：

2018宁夏网络赛H题：
题意：
给你类似剪刀石头布游戏的五种手势和十种克制关系。每种手势克制其他两种手势并被另外两种手势克制。
给你两个字符串分别表示A，B的手势顺序，A长B短，你可以随意挪动B相对于A的位置，求B最多能赢多少次？
解题思路：我们固定A串，发现B串只有固定次数的匹配方式，而如果我们把短串逆置，我们就可以用一定的指数来表示两串的某种位置关系。
例如
原串 abcdef
短串 ghi
将短串逆置
原串 a b c d e f
逆置串 i h g
我们给他们标上次幂
$a*x^{5}quad quad b*x^{4}quad quad c*x^{3}quad quad d*x^{2}quad quad e*xquad quad f$
$i*x^{2}quad quad h*xquad quad quad g$
如果将两个多项式相乘，那么对于 $x^{5}$ 则是短串与原串前三个位置匹配的结果。
$x^{4}$ 则是短串与原串2-4匹配的结果。依此类推…
那么我们如何让多项式乘法的系数来表示赢多少次呢，仔细考虑一下我们会发现我们不能一次算出所有赢的情况，但是我们可以把所有种赢的情况分别进行多项式乘法，最后取和。
如果a>b&&a>c
那么我们将短串a标记为1，长串的b,c标记为1，其余都标记为0。
我们会发现只有a与b,c进行匹配的乘法才会有答案，所以进行一次FFT之后，各个指数的系数即为该匹配位置中，a赢b,c的次数，同理，我们可以算出b中所有赢的方式赢得次数，最后对所有匹配位置赢得次数之和取max，即为本题最终解。
附上代码：

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

const double pi=acos(-1.0);
struct complex
{
    double x,y;
    complex(double X=0,double Y=0)
    {
        x=X,y=Y;
    }
}a[maxn],b[maxn];
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);}
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(complex a[maxn],int opt)
{
    for (int i=0;i<n;++i)
        if (i<R[i]) swap(a[i],a[R[i]]);
    for (int k=1;k<n;k<<=1)
    {
        complex wn=complex(cos(pi/k),opt*sin(pi/k));
        for (int i=0;i<n;i+=(k<<1))
        {
            complex w=complex(1,0);
            for (int j=0;j<k;++j,w=w*wn)
            {
                complex x=a[i+j],y=w*a[i+j+k];
                a[i+j]=x+y,a[i+j+k]=x-y;
            }
        }
    }
}
void calc(int opt)
{
    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].x/n+0.5)*opt;//对于每种匹配位置累计赢的次数
    }
}
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;
    for (int i=0;i<n;++i)
        R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
    //S>P&&S>L
    for (int i=0;i<T;++i) f[i]=(t[i]=='P'||t[i]=='L')?1:0;
    for (int i=0;i<S;++i) g[i]=(s[i]=='S')?1:0;
    for (int i=0;i<=n;++i) a[i]=complex(0,0),b[i]=complex(0,0);
    for (int i=0;i<T;++i) a[i].x=f[i];
    for (int i=0;i<S;++i) b[i].x=g[i];
    calc(1);
    //P>R&&P>K
    for (int i=0;i<T;++i) f[i]=(t[i]=='K'||t[i]=='R')?1:0;
    for (int i=0;i<S;++i) g[i]=(s[i]=='P')?1:0;
    for (int i=0;i<=n;++i) a[i]=complex(0,0),b[i]=complex(0,0);
    for (int i=0;i<T;++i) a[i].x=f[i];
    for (int i=0;i<S;++i) b[i].x=g[i];
    calc(1);
    //R>S&&R>L
    for (int i=0;i<T;++i) f[i]=(t[i]=='S'||t[i]=='L')?1:0;
    for (int i=0;i<S;++i) g[i]=(s[i]=='R')?1:0;
    for (int i=0;i<=n;++i) a[i]=complex(0,0),b[i]=complex(0,0);
    for (int i=0;i<T;++i) a[i].x=f[i];
    for (int i=0;i<S;++i) b[i].x=g[i];
    calc(1);
    //L>P&&L>K
    for (int i=0;i<T;++i) f[i]=(t[i]=='P'||t[i]=='K')?1:0;
    for (int i=0;i<S;++i) g[i]=(s[i]=='L')?1:0;
    for (int i=0;i<=n;++i) a[i]=complex(0,0),b[i]=complex(0,0);
    for (int i=0;i<T;++i) a[i].x=f[i];
    for (int i=0;i<S;++i) b[i].x=g[i];
    calc(1);
    //K>R&&K>S
    for (int i=0;i<T;++i) f[i]=(t[i]=='R'||t[i]=='S')?1:0;
    for (int i=0;i<S;++i) g[i]=(s[i]=='K')?1:0;
    for (int i=0;i<=n;++i) a[i]=complex(0,0),b[i]=complex(0,0);
    for (int i=0;i<T;++i) a[i].x=f[i];
    for (int i=0;i<S;++i) b[i].x=g[i];
    calc(1);

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

FFT第二题：

BZOJ4259(感谢M哥的权限号
题意：
给你两个字符串原串和匹配串，串里只包含小写字母和 $*$ , $*$ 可以表示任意字符，问匹配串在原串不同位置中出现多少次，起始位置不同即不同位置。
解题思路：
我们还是利用上一题逆置短串构造卷积的方法，现在要思考的是怎么让系数表示该位置是否可以匹配呢？
仔细思考一下，两个位置匹配无非三种情况，对于所有的位置下标
1. $str1[i]=str2[i]$
2. $str1[i]=$ $*$
3. $str2[i]=$ $*$
以上三种情况满足任意一项即可，
所以我们可以用 $(str1[i]-str2[i])^{2}*str1[i]*str2[i]$ 来表示两个字符串某位置的匹配情况，如果某项系数为0，则代表位置匹配。至于第一项为什么要平方，是因为我们保证所有项都是正数，以免相加结果为0影响答案。
现在我们展开上面的式子。
$quad (str1[i]-str2[i])^{2}*str1[i]*str2[i]$
$=((str1[i])^{3}*str2[i]-2*(str1[i])^{2}*(str2[i])^{2}+str1[i]*(str2[i])^{3})$
现在我们如果反转str2,就可以得到卷积
$F(i)=((str1[i])^{3}*str2[n-i]-2*(str1[i])^{2}*(str2[n-i])^{2}+str1[i]*(str2[n-i])^{3})$
然后分三段进行FFT再累计结果就可以。
附上代码：

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

const double pi=acos(-1.0);
struct complex
{
    double x,y;
    complex(double X=0,double Y=0)
    {
        x=X,y=Y;
    }
}a[maxn],b[maxn];
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);}
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(complex a[maxn],int opt)
{
    for (int i=0;i<n;++i)
        if (i<R[i]) swap(a[i],a[R[i]]);
    for (int k=1;k<n;k<<=1)
    {
        complex wn=complex(cos(pi/k),opt*sin(pi/k));
        for (int i=0;i<n;i+=(k<<1))
        {
            complex w=complex(1,0);
            for (int j=0;j<k;++j,w=w*wn)
            {
                complex x=a[i+j],y=w*a[i+j+k];
                a[i+j]=x+y,a[i+j+k]=x-y;
            }
        }
    }
}
void calc(int opt)
{
    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].x/n+0.5)*opt;
}
int main()
{
    scanf("%d%d",&S,&T);
    scanf("%s%s",s,t);
    for (int i=0;i<S/2;++i) swap(s[i],s[S-i-1]);

    for (int i=0;i<T;++i) f[i]=(t[i]=='*')?0:(t[i]-'a'+1.0);//预处理*的位置
    for (int i=0;i<S;++i) g[i]=(s[i]=='*')?0:(s[i]-'a'+1.0);

    m=S+T-2;
    for (n=1;n<=m;n<<=1) ++L;
    for (int i=0;i<n;++i)
        R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));

    for (int i=0;i<=n;++i) a[i]=complex(0,0),b[i]=complex(0,0);
    for (int i=0;i<T;++i) a[i].x=f[i]*f[i]*f[i];//按照卷积表达式进行修改
    for (int i=0;i<S;++i) b[i].x=g[i];
    calc(1);//对于每一项进行FFT，参数即为卷积表达式中的系数
    for (int i=0;i<=n;++i) a[i]=complex(0,0),b[i]=complex(0,0);
    for (int i=0;i<T;++i) a[i].x=f[i]*f[i];
    for (int i=0;i<S;++i) b[i].x=g[i]*g[i];
    calc(-2);
    for (int i=0;i<=n;++i) a[i]=complex(0,0),b[i]=complex(0,0);
    for (int i=0;i<T;++i) a[i].x=f[i];
    for (int i=0;i<S;++i) b[i].x=g[i]*g[i]*g[i];
    calc(1);

    for (int i=S-1;i<T;++i)
        if (!F[i]) ans[++ans[0]]=i-S+2;//每种指数的起始位置
    printf("%dn",ans[0]);
    for (int i=1;i<=ans[0];++i) printf("%d%c",ans[i]," n"[i==ans[0]]);
}

FFT第三题：

2016香港网络赛
题意：
给你一个大小为2e5的数组，问有多少对a[i]+a[j]==a[k],i,j,k两两不同。
解题思路：
我们可以把这题换一种题意理解，去数组中所有的i,j(i!=j)a[i],令a[i],a[j]相加，问有多少个答案在数组中。所以就变成了第二类问题，统计两个数组所有加法结果的出现次数，而本题则是数组与本身的所有加法。首先考虑，把数组中所有元素转化为指数，该元素的出现次数转换为系数，例如
原数组为 1 1 1 2 3 3 4 4
则可看成 $3*x^{1}+x^{2}+2*x^{3}+2*x^{4}$
而如果两个相同的多项式相乘，得到的指数即为所有可能加法的结果，而对应的系数即为加法结果出现次数。这个题大概就解决了，但是这题细节比较多，所以我们要注意以下几点。
1.原题中会有负数出现，为避免对结果产生影响，我们将所有数均加上一个相同的大数，以保证所有指数都为正数。
2.统计结果中会有某元素与本身相加的结果，所以对每个元素x的结果要减去2*x两次。
3.设我们加的值为LOW，则统计式可改为 $a+LOW+b+LOW==c+2*LOW$ ,所以我们要统计每个值 $+2*LOW$ 出现的次数。
4.考虑0对于其他元素的影响，
对于不等于0的元素x来说
$0+x=x,x+0=x$
所以要减去0的个数的二倍。
对于等于0的元素0来说
$0+0=0,0+0=0$
但是有本元素在左在右两个位置的结果是要统计的，所以只要减去0的个数减一的二倍。

对FFT的结果用以上四种条件进行约束，最终结果即为本题答案。
附上代码：

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
using namespace std;
#define maxn 800005
typedef long long ll;
const int LOW=50000;
const double pi=acos(-1.0);
struct complex
{
    double x,y;
    complex(double X=0,double Y=0)
    {
        x=X,y=Y;
    }
}a[maxn],b[maxn];
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);}
int S,T,n,m,L,R[maxn],zero;
long long F[maxn],num[maxn],cnt[maxn];
void FFT(complex a[maxn],int opt)
{
    for (int i=0;i<n;++i)
        if (i<R[i]) swap(a[i],a[R[i]]);
    for (int k=1;k<n;k<<=1)
    {
        complex wn=complex(cos(pi/k),opt*sin(pi/k));
        for (int i=0;i<n;i+=(k<<1))
        {
            complex w=complex(1,0);
            for (int j=0;j<k;++j,w=w*wn)
            {
                complex x=a[i+j],y=w*a[i+j+k];
                a[i+j]=x+y,a[i+j+k]=x-y;
            }
        }
    }
}
void calc(int opt)
{
    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=0;i<n;++i)//对于统计两数组加法的都要完全统计0-n
    {
         F[i]+=(long long)(a[i].x/n+0.5)*opt;
    }
}
int main()
{
    scanf("%d",&S);
    T=S;
    for(int i=0;i<S;i++)
    {
        scanf("%lld",&num[i]);
        if(num[i]==0) zero++;
        cnt[num[i]+LOW]++;//离散化
    }
    for (n=1;n<=200000;n<<=1) ++L;
    for (int i=0;i<n;++i)
        R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
    for (int i=0;i<=n;++i) a[i]=complex(1.0*cnt[i],0.0),b[i]=complex(1.0*cnt[i],0.0);
    calc(1);//传递的参数为多项式乘法的常数可以为正负

    for(int i=0;i<T;i++)//去除本身与本身相加的结果
    {
        F[(num[i]+LOW)*2]--;
    }
    ll pp = 0;
    for(int i=0;i<T;i++)
    {
        pp+=F[num[i]+LOW*2];
        if(num[i]==0)
            pp-=(zero-1)*2;//元素为0的情况
        else
            pp-=zero*2;//元素非0的情况
    }
    printf("%lldn",pp);
    return 0;
}