FFT:从入门到沉迷

终于学会了$FFT$并无法自拔

啥是FFT

$FFT$指快速傅里叶变换

在$OI$中的应用是$O(nlogn)$计算函数卷积

人话：多项式乘法

多项式毒瘤模板题的万恶之源

系数表达式与点值表达式

系数表达式就是平常的表示方法

$$f(x)=\sum _{i=0}^n{a}_i{x}^i$$

点值表达式就是取函数上$n$个不同的点来表示

如：多项式$f(x)=x^3-2x^2-x+1$

点值表达式可以为：$(-1,-1)(0,1)(1,-1)$

由代数基本定理，两者等价

初步推导

为什么要引入点值表达式？

如果相乘的多项式选取了相同的$x$，就可以$O(n)$计算乘积

现在的问题是：由系数转点值，乘完后再转系数

然而系数转点值需要取$n$个$x$,然后求值，复杂度$O(n^2)$

更尴尬的是点值转系数高斯消元$O(n^3)$

看来要在取数上做手脚了

复平面与单位根

定义$i^2=-1$，则$a+bi$表示一个复数

不需要纠结$i$是个啥，只需要保留，遇到$i^2$就化成$-1$就行了

定义复平面上$(a,b)$表示$a+bi$,也可看成向量

有个奇妙的性质：复数相乘，相当于转角（$x$正半轴逆时针旋转）相加，长度相乘

定义$n$的单位根$w_n$满足$w_n^n=1$

这样的单位根有$n$个，为$w_n^0,w_n^1,w_n^2,...,w_n^{n-1}$

他们将平分整个复平面

折半引理

$w_n^k=w_{2n}^{2k}$

之所以叫折半而不叫倍增，主要是分数太难看所以转了一下

画个图就能证明

快速傅里叶变换

为表述方便，我们把不足的$n$增加次数凑成$2$的整数次幂。注意$n$为长度，最高次为$n-1$

脑子一热，不如我们把$n$个单位根作$x$吧！

开始推式子了

$$f(x)=\sum _{i=0}^{n-1}{a}_i{x}^i$$

闲着没事按奇偶分类

$$f(x)=\sum _{i=0}^{\frac{n}{2}}{a}_{2i}{x}^{2i}+\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i+1}{x}^{2i+1}$$

太丑了

定义

$$f_1(x)=\sum _{i=0}^{\frac{n}{2}}{a}_{2i}{x}^i$$

$$f_2(x)=\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i+1}{x}^i$$

那么有

$$f(x)=f_1(x^2)+x{f}_2(x^2)$$

对于$k\le \frac{n}{2}$,代入$x=w_n^k$

$$f(w_n^k)=f_1(w_n^{2k})+w_n^k{f}_2(w_n^{2k})$$

同时

$$f(w_n^{k+\frac{n}{2}})=f_1(w_n^{2k+n})+w_n^{k+\frac{n}{2}}{f}_2(w_n^{2k+n})$$

显然$w_n^n=1,w_{2n}^n=-1$

所以

$$f(w_n^{k+\frac{n}{2}})=f_1(w_n^{2k})-w_n^k{f}_2(w_n^{2k})$$

我们发现两个式子结构是相同的

根据折半引理

$$f_1(w_n^{2k})=f_1(w_{\frac{n}{2}}^k)$$

$$f_2(w_n^{2k})=f_2(w_{\frac{n}{2}}^k)$$

故

$$f(w_n^k)=f_1(w_{\frac{n}{2}}^k)+w_n^k{f}_2(w_{\frac{n}{2}}^k)$$

$$f(w_n^{k+\frac{n}{2}})=f_1(w_{\frac{n}{2}}^k)-w_n^k{f}_2(w_{\frac{n}{2}}^k)$$

我们发现规模不断减半，然后可以愉快的递归了。

复杂度$O(NlogN)$

逆变换

现在我们转成了点值，进行了乘法，该转回系数了。

假设原多项式系数是$a$，将$N$个单位根代入得到了$N$个复数$b$

即：$b_k={\sum }_{i=0}^{n-1}{a}_i{w}_n^{ki}$

现在我们要求$a$

把$b$再插回去试试

我们倒着插单位根

即$c_k$为多项式$$\sum _{i=0}^{n-1}{b}_i{w}_n^{-ki}$$

代入$b$

$$c_k=\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}{a}_j{w}_n^{ij}{w}_n^{-ki}$$

把$w_n$合并

$$c_k=\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}{a}_j{w}_n^{ij-ki}$$

$$=\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}{a}_j{w}_n^{i(j-k)}$$

$a_j$可以放前面来

$$c_k=\sum _{j=0}^{n-1}{a}_j\sum _{i=0}^{n-1}{w}_n^{i(j-k)}$$

考虑$$S=\sum _{i=0}^{n-1}{w}_n^{i(j-k)}$$

当$j=k$,显然$S=n$

当$j̸=k$

$$w_n^{j-k}S=\sum _{i=1}^n{w}_n^{i(j-k)}$$

$$(w_n^{j-k}-1)S=w_n^{n(j-k)}-w_n^{j-k}=w_n^{j-k}-w_n^{j-k}=0$$

所以$S=0$

这说明只有$j=k$时对答案有贡献

此时

$$c_k=a_{k}n$$

$$a_k=\frac{c_k}{n}$$

$FFT$至此就结束了。

代码：

#include <iostream>
#include <cstdio>
#include <cmath>
#include <cctype>
#define MAXN 2000005
const double PI=acos(-1.0);
using namespace std;
inline int read()
{
    int ans=0,f=1;
    char c=getchar();
    while (!isdigit(c))
    {
        if (c=='-')
            f=-1;
        c=getchar();
    }
    while (isdigit(c))
    {
        ans=(ans<<3)+(ans<<1)+(c^48);
        c=getchar();
    }
    return f*ans;
}
struct complex
{
    double x,y;
    complex(double x=0,double y=0):x(x),y(y){}
}a[MAXN],b[MAXN];
complex operator +(const complex& a,const complex& b){return complex(a.x+b.x,a.y+b.y);}
complex operator -(const complex& a,const complex& b){return complex(a.x-b.x,a.y-b.y);}
complex operator *(const complex& a,const complex& b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
void fft(int lim,complex *a,int type)
{
    if (lim==1)
        return;
    complex a1[lim>>1],a2[lim>>1];
    for (int i=0;i<=lim;i+=2)
        a1[i>>1]=a[i],a2[i>>1]=a[i+1];
    fft(lim>>1,a1,type);
    fft(lim>>1,a2,type);
    complex wn(cos(2.0*PI/lim),type*sin(2.0*PI/lim)),w(1,0);
    for (int i=0;i<(lim>>1);i++,w=w*wn)
    {
        a[i]=a1[i]+w*a2[i];
        a[i+(lim>>1)]=a1[i]-w*a2[i];
    }
}
int main()
{
    int n,m;
    n=read(),m=read();
    for (int i=0;i<=n;i++)//注意下标 
        a[i].x=read();
    for (int i=0;i<=m;i++)
        b[i].x=read();
    int lim=1;
    while (lim<=n+m)
        lim<<=1;
    fft(lim,a,1);
    fft(lim,b,1);
    for (int i=0;i<=lim;i++)	
        a[i]=a[i]*b[i];
    fft(lim,a,-1);
    for (int i=0;i<=n+m;i++)
        printf("%d ",(int)(a[i].x/lim+0.5));
    return 0;
}

迭代实现

什么？你没过？

观察一下原来的代码 我们发现数组开小了

在递归函数中，每一层都会开一个数组，这造成效率极低。

我们观察递归的全过程

（二进制）

000 001 010 011 100 101 110 111
000 010 100 110|001 011 101 111
000 100|010 110|001 101|011|111

我们发现，递归完后实际上是将二进制翻转

那我们可以模拟递归的过程，算出最后的下标，再一层一层推回去

如何计算翻转后的数？

①将第二位及以上的翻转并右移1位

②将第一位移到最高位

可以$O(N)$求出

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cctype>
#include <cmath>
#define MAXN 5000005
const double Pi=acos(-1.0);
using namespace std;
inline int read()
{
    int ans=0,f=1;
    char c=getchar();
    while (!isdigit(c)){if (c=='-') f=-1;c=getchar();}
    while (isdigit(c)) ans=(ans<<3)+(ans<<1)+(c^48),c=getchar();
    return f*ans;
}
struct complex{double x,y;complex(const double& x=0,const double& y=0):x(x),y(y){}}a[MAXN],b[MAXN];
inline complex operator +(const complex& a,const complex& b){return complex(a.x+b.x,a.y+b.y);}
inline complex operator -(const complex& a,const complex& b){return complex(a.x-b.x,a.y-b.y);}
inline complex operator *(const complex& a,const complex& b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
int lim=1,l,r[MAXN];
void fft(complex* a,int type)
{
    for (int i=0;i<lim;i++) if (i<r[i]) swap(a[i],a[r[i]]);// 计算最后的序列。<是为了避免换回来，写>也一样
    for (int len=2;len<=lim;len<<=1)//合并出的长度
    {
        int mid=(len>>1);//被合并的长度
        complex Wn(cos(Pi/mid),type*sin(Pi/mid));//单位根
        for (int s=0;s<=lim;s+=len)//枚举每一段的长度
        {
            complex w(1,0);
            for (int k=0;k<mid;k++,w=w*Wn)
            {
                complex x=a[s+k],y=w*a[s+mid+k];//一个小优化，避免大量复数乘法运算
                a[s+k]=x+y,a[s+mid+k]=x-y;
            }
        }
    }
}
int main()
{
    int n,m;
    n=read(),m=read();
    for (int i=0;i<=n;i++) a[i].x=read();
    for (int i=0;i<=m;i++) b[i].x=read();
    while (lim<=n+m) lim<<=1,++l;
    for (int i=0;i<lim;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));//r[i]表示i翻转后的数
    fft(a,1),fft(b,1);
    for (int i=0;i<=lim;i++) a[i]=a[i]*b[i];
    fft(a,-1);
    for (int i=0;i<=n+m;i++) printf("%d ",(int)(a[i].x/lim+0.5));
    return 0;
}

凡是两个数组一一对应然后进行各种运算的，都可以考虑$FFT$

甚至可以用来匹配字符串