FFT

前置知识

复数的指数形式

\(a+bi=re^{i\theta} (r是长度，\theta 是在复平面跟x轴的角度), r=\sqrt{a^2+b^2},tan(\theta)=\frac{b}{a}\)

单位根

\(x^n=1\)在复数域的根称为n次单位根，n次单位根有n个，形式为\((\omega_n)^k=\omega^k_n=e^{i\frac{2k\pi}{n}}=(cos(\frac{2k\pi}{n}),sin(\frac{2k\pi}{n})i)\)

单位根的性质

\(\omega^k_n=\omega^{2k}_{2n}\)
\(\omega^{k+n}_{2n}=\omega_{2n}^k*\omega_{2n}^n=-\omega^{k}_{2n}\)

多项式的系数表示和点值表示

假设\(f(x)\)为一个n次多项式，则\(f(x)\)的系数表示为\(f(x)=a_0+a_1x+a_2x^2+···+a_{n}x^{n}\)
\(f(x)\)的点值表示为\((x_n,f(x_n)),(x_{n-1},f(x_{n-1})),···,(x_0,f(x_0))\)

• n+1个点可以表示一个n次多项式
• 在点值表示下n次多项式的乘法复杂度为\(O(n)\)
  假设\(h(x)=f(x)g(x),并且知道(x_i,f(x_i)),(x_i,g(x_i)), 1\leq i\leq n\)
  则知道了\((x_i,h(x_i))=(x_i,f(x_i)g(x_i))\)
  就\(O(n)\)知道了h(x)的点值表示

那么对于多项式乘法

F(x)和G(x)的系数表示形式--->F(x),G(x)的点值表示形式---\(O(n)\)--->H(x)的点值形式--->H(x)的系数形式
关键在于系数形式和点值形式的相互转化，FFT就充当这个角色，在\(O(nlogn)\)就能相互转化

FFT(快速傅里叶变换)

DFT(离散傅里叶变换)

将多项式\(A(x)=a_0+a_1x+a_2x^2+···+a_{n-1}x^{n-1}\)转化为点值形式\((\omega_n^k,A(\omega_n^k)),(k=0,1,2,···,n-1)\)
假设\(n=2^l\)
\(A(x)=(a_0+a_2x^2+···+a_{n-2}x^{n-2})+(a_1x+a_3x^3+···+a_{n-1}x^{n-1})\)
令
\(B(x)=a_0+a_2x+···+a_{n-2}x^{\frac{n}{2}-1}\)
\(C(x)=a_1+a_3x+···+a_{n-1}x^{\frac{n}{2}-1}\)
\(A(x)=B(x^2)+xC(x^2)\)

对于\(0\leq k<\frac{n}{2}\)

\(A(\omega_n^k)=B(\omega_n^{2k})+\omega_n^kC(\omega_n^{2k})\)
\(\:\:\:\:\:\:\:\:\:\:\:\:=B(\omega_{\frac{n}{2}}^k)+\omega_n^kC(\omega_{\frac{n}{2}}^k)\)

对于\(\frac{n}{2}\leq k\)
\(A(\omega_n^{k+\frac{n}{2}})=B(\omega_n^{2k+n})+\omega_n^{\frac{n}{2}+k}C(\omega_n^{2k+n})\)
\(\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:=B(\omega_n^{2k})+\omega_n^{\frac{n}{2}+k}C(\omega_n^{2k})\)

\(\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:=B(\omega_{\frac{n}{2}}^k)-\omega_n^{k}C(\omega_{\frac{n}{2}}^{k})\)

这样就得到了一个分治递归做法，但常数大
蝴蝶变换非递归做法常数小
简单来说，分治做法是由上到下分治到最底层再到上合并的递归过程
蝴蝶变换直接通过二进制就能得到最底层的分治顺序，然后再向上合并，常数小

蝴蝶变换

会发现序列的初始顺序的每个数二进制与分治最终的顺序的每个数的二进制是相反的
具体实现：
一个数的二进制反转为这个数除2的反转除2（得到后（l-1）位），如果这个数第一位为1，那么把最后一位加1。

点击查看代码

void get(int bit){
    for(int i=0;i<(1<<bit);i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    return ;
}

###IDFT(逆离散傅里叶变换) 将多项式的点值表示$(\omega_n^k,b_k),(k=0,1,2,···,n-1)$转化为其系数表示$A(x)=a_0+a_1x+a_2x^2+···+a_{n-1}x^{n-1}$ 将点值$(\omega_n^k,b_k)$代入$A(x)$,则可以转化为矩阵相乘 $$\begin{pmatrix} \omega_n^{0\cdot 0}&\omega_n^{0\cdot 1}&···&\omega_n^{0\cdot (n-1)}\\ \omega_n^{1\cdot 0}&\omega_n^{1\cdot1}&···&\omega_n^{1\cdot (n-1)}\\ ·&·&&·\\ ·&·&&·\\ \omega_n^{(n-1)\cdot 0}&\omega_n^{(n-1)\cdot1 }&···&\omega_n^{(n-1)\cdot(n-1)}\\ \end{pmatrix} \cdot\begin{pmatrix} a_0\\a_1\\·\\·\\a_{n-1} \end{pmatrix}=\begin{pmatrix} b_0\\b_1\\·\\·\\b_{n-1} \end{pmatrix}$$ 设最左边的矩阵为$\Omega$,它的逆元为$\Omega^{-1}$ 则 $$\begin{pmatrix} a_0\\a_1\\·\\·\\a_{n-1} \end{pmatrix}=\Omega^{-1}\cdot\begin{pmatrix} b_0\\b_1\\·\\·\\b_{n-1} \end{pmatrix}$$

则IDFT就是求出\(\Omega^{-1}\)
设

\[\overline{\Omega}=\begin{pmatrix} \omega_n^{-0\cdot 0}&\omega_n^{-0\cdot 1}&···&\omega_n^{-0\cdot (n-1)}\\ \omega_n^{-1\cdot 0}&\omega_n^{-1\cdot1}&···&\omega_n^{-1\cdot (n-1)}\\ ·&·&&·\\ ·&·&&·\\ \omega_n^{-(n-1)\cdot 0}&\omega_n^{-(n-1)\cdot1 }&···&\omega_n^{-(n-1)\cdot(n-1)}\\ \end{pmatrix}\]

则令\(M=\overline{\Omega}\cdot\Omega\)
\(M_{ij}=\omega_n^{-i\cdot 0}\omega_n^{0\cdot j}+\omega_n^{-i\cdot 1}\omega_n^{1\cdot j}+···+\omega_n^{-i\cdot (n-1)}\omega_n^{(n-1)\cdot j}\)
\(\:\:\:\:\:\:\:\:=\omega_n^{(j-i)\cdot0}+ \omega_n^{(j-i)\cdot1}+···+\omega_n^{(j-i)\cdot(n-1)}\)
所以，

\[ M_{ij}=\left\{ \begin{aligned} n & , & j-i=0 \: , \:\omega_n^0=1 \\ \frac{1-(\omega_n^{j-i})^n}{1-\omega_n^{j-i}}=0（(\omega_n^{j-i})^n=1） &, & j-i\not=0\:,\:\omega_n^{i-j}\not=1 \end{aligned} \right. \]

则

\[M=\begin{pmatrix} n&0&0&····&0\\ 0&n&0&····&0\\ ·&·&·&&·\\ ·&·&·&&·\\ 0&0&0&····&n\\ \end{pmatrix}=nI(I为单位矩阵)\]

因此\(\overline{\Omega}\cdot\Omega=nI\)
\(\Omega^{-1}=\frac{1}{n}\overline{\Omega}\)
由此我们可以得到

\[\begin{pmatrix} a_0\\a_1\\·\\·\\a_{n-1} \end{pmatrix}=\frac{1}{n}\overline{\Omega}\cdot\begin{pmatrix} b_0\\b_1\\·\\·\\b_{n-1} \end{pmatrix}=\frac{1}{n}\begin{pmatrix} b_0\omega_n^{-0\cdot0}+b_1\omega_n^{-0\cdot1}+···+b_{n-1}\omega_n^{-0\cdot(n-1)}=B(\omega_n^{-0})\\ b_0\omega_n^{-1\cdot0}+b_1\omega_n^{-1\cdot1}+···+b_{n-1}\omega_n^{-1\cdot(n-1)}=B(\omega_n^{-1})\\ ·\\·\\ b_0\omega_n^{-(n-1)\cdot0}+b_1\omega_n^{-(n-1)\cdot1}+···+b_{n-1}\omega_n^{-(n-1)\cdot(n-1)}=B(\omega_n^{-(n-1)})\\ \end{pmatrix}\]

相当于给定\(B(x)=b_0+b_1x+···+b_{n-1}x^{n-1}\),求点值\(B(\omega_n^{-k}),(k=0,1,2,···,n-1)\),这样就可以用DFT来实现

总结

FFT就是对于\(H(x)=F(x)G(x)\),先分别做一遍DFT求出\(F(x)和G(x)的点值表示f(x),g(x)\),\(H(x)的点值表示h(x)=f(x)g(x)\),在对于\(h(x)\)做一遍DFT求出\(H(x)\)的系数

模板

点击查看代码

#include<algorithm>
#include<iostream>
#include<complex>
#include<cstdlib>
#include<cstring>
#include<cstdio>
#include<cmath>
#include<queue>
#define ll long long 
using namespace std;
const int maxn=8000000+10101;
const double pi=acos(-1); //圆周率  
inline int read(){
    int x=0,f=1;char ch=getchar();
    for(;!isdigit(ch);ch=getchar())if(ch=='-')f=-1;
    for(;isdigit(ch);ch=getchar())x=(x<<3)+(x<<1)+ch-'0';
    return x*f; 
}
typedef complex<double> cd; //复数类的定义 
int rev[maxn];              //rev[i] 表示 i 二进制反转的 结果，也就是第i个系数在进行奇偶次项分治后的位置 
void get(int bit){//二进制翻转 
    for(int i=0;i<(1<<bit);i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    return ;
}
int n,m;
cd a[maxn],b[maxn];
void dft(cd *u,int val){
    for(int i=0;i<n;i++)if(i<rev[i])swap(u[i],u[rev[i]]);
    //保证每个数只被交换一次。设x=1(001),则rev[x]=4(100),当从前往后进行交换，
    //第一次扫到x时交换x,rev[x],x就到就到后面了，当再次扫到x再次交换，
    //x,rev[x]的顺序与原始状态无改变 
    for(int i=1;i<n;i<<=1){//枚举待合并区间的一半 
        cd wn(cos(pi/i),val*sin(pi/i));
        //计算单位根，因为i是待合并区间的一半，所以2*pi/(i*2)》pi/i 
        for(int j=0;j<n;j+=(i<<1)){//对于每一区间 
            cd w(1,0);//每一块都是一个独立序列,都是以零次方位为起始的,可以当做幂 
            for(int k=0;k<i;k++,w*=wn){//蝴蝶操作处理这一区间 
                cd x=u[j+k],y=w*u[j+k+i];
                u[j+k]=x+y;
                u[j+k+i]=x-y;//利用单位根的性质，O(1)求出另一半 
            }
        }
    }
    return ;
}
int main(){
    n=read();m=read();
    for(int i=0;i<=n;i++)a[i]=read();
    for(int i=0;i<=m;i++)b[i]=read();
    int len=0;m+=n;
    for(n=1;n<=m;n<<=1)len++;
    //找到第一个二的整数次幂使得其可以容纳这两个数的乘积 
    //因为在蝴蝶操作过程一定要求每层都可以分成两大小相等的部分，所以多项式最高次项一定是
    get(len);dft(a,1);dft(b,1);
    for(int i=0;i<=n;i++)a[i]*=b[i];
    dft(a,-1);//IDFT
    for(int i=0;i<=m;i++)printf("%d ",(int)(a[i].real()/n+0.5)); //注意精度误差 
    return 0;   
}

应用

1.A * B Problem Plu

点击查看代码

#include<functional>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<complex>
#include<string>
#include<cstdio>
#include<vector>
#include<cmath>
#include<queue>
#include<deque>
#define ll long long 
using namespace std;
const int maxn=2e5+10101;
const int MOD=1e9+7;
const int inf=2147483647;
const double pi=acos(-1);
int read(){
    int x=0,f=1;char ch=getchar();
    for(;!isdigit(ch);ch=getchar())if(ch=='-')f=-1;
    for(;isdigit(ch);ch=getchar())x=x*10+ch-'0';
    return x*f;
}
typedef complex<double> cd;
int n,m,rev[maxn],len,c[maxn],d[maxn];
cd a[maxn],b[maxn];
void get(int bit){
    for(int i=0;i<(1<<bit);i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
}
void dft(cd *u,int val){
    for(int i=0;i<n;i++)if(i<rev[i])swap(u[i],u[rev[i]]);
    for(int i=1;i<n;i<<=1){
        cd wn(cos(pi/i),val*sin(pi/i));
        for(int j=0;j<n;j+=(i<<1)){
            cd w(1,0);
            for(int k=0;k<i;k++,w*=wn){
                cd x=u[j+k],y=w*u[j+k+i];
                u[j+k]=x+y;
                u[j+k+i]=x-y;
            }
        }
    }
    return ;
}
int main(){
    char ch[maxn];cin>>ch;n=strlen(ch);
    for(int i=0;i<n;i++)a[i]=(cd)(ch[n-i-1]-'0');
    cin>>ch;m=strlen(ch);
    for(int i=0;i<m;i++)b[i]=(cd)(ch[m-i-1]-'0');
    m+=n;for(n=1;n<=m;n<<=1)len++;
    get(len);dft(a,1);dft(b,1); 
    for(int i=0;i<(1<<len);i++)a[i]*=b[i];
    dft(a,-1);
    for(int i=0;i<m;i++)c[i]=(int)(a[i].real()/n+0.5);
    int jin=0,tot=0;
    for(int i=0;i<m;i++){
        d[i]=(c[i]+jin)%10;
        if(c[i]+jin>=10)jin=(c[i]+jin)/10;
        else jin=0;
    }
    bool fa=false;
    for(int i=m-1;i>=0;i--){
        if(d[i]!=0)fa=true;
        if(fa)printf("%d",d[i]);
    }
    return 0;
}

2.TSUM - Triple Sums 题解

3.Fuzzy Search 题解