快速傅立叶变换

前置知识

离散傅立叶变换（DFT）

快速傅立叶变换（FFT）

$X(k)={\sum }_{n=0}^{N-1}x(n)e^{-\frac{2\pi i}{N}nk}$
$x(k)=\frac{1}{N}{\sum }_{k=0}^{N-1}x(n)e^{\frac{2\pi i}{N}nk}$
FFT简单来说就是加速DFT的
显然计算一次DFT的时间复杂度是$O(N^2)$
而FFT可以让他降为$O(N\log N)$

推导

记$W_N=e^{-\frac{2\pi i}{N}}$，也叫旋转因子
则$X(k)={\sum }_{n=0}^{N-1}x(k)W_N^{nk}$
然后我们推导一下$W_N^{nk}$的性质
1.$W_N^{rN}=1(r=0,1,\cdots )$
2.周期性
  $W_N^{n(k+rN)}=W_N^{nk}$
3.对称性
  $W_N^{nk+\frac{N}{2}}=-W_N^{nk}$
4.$W_{2N}^{2nk}=W_N^{nk}$

接下来我们正式推导
假设$N$是$2$的整数幂
$$\begin{array}{rl}& X(k)\\ & =\sum _{n=0}^{N-1}x(k)W_N^{nk}\\ & =\sum _{r=0}^{\frac{N}{2}-1}x(2r)W_N^{2rk}+\sum _{r=0}^{\frac{N}{2}-1}x(2r+1)W_N^{k(2r+1)}\\ & =\sum _{r=0}^{\frac{N}{2}-1}x(2r)W_N^{2rk}+W_N^k\sum _{r=0}^{\frac{N}{2}-1}x(2r+1)W_N^{2rk}\\ & =\sum _{r=0}^{\frac{N}{2}-1}x(2r)W_{\frac{N}{2}}^{rk}+W_N^k\sum _{r=0}^{\frac{N}{2}-1}x(2r+1)W_{\frac{N}{2}}^{rk}\end{array}$$
记
  $G(k)={\sum }_{r=0}^{\frac{N}{2}-1}x(2r)W_{\frac{N}{2}}^{rk}$
  $H(k)={\sum }_{r=0}^{\frac{N}{2}-1}x(2r+1)W_{\frac{N}{2}}^{rk}$
则$X(k)=G(k)+W_N^{k}H(k)$
利用周期性
$$G(k+\frac{N}{2})=\sum _{r=0}^{\frac{N}{2}-1}x(2r)W_{\frac{N}{2}}^{r(k+\frac{N}{2})}=G(k)$$
同理
$$H(k+\frac{N}{2})=H(k)$$
于是我们得到
$$X(k+\frac{N}{2})=G(k+\frac{N}{2})+W_N^{k+\frac{N}{2}}H(k+\frac{N}{2})=G(k)-W_N^{k}H(k)$$
对比
$$X(k)=G(k)+W_N^{k}H(k)$$
我们发现，间隔$\frac{N}{2}$的$X$，只差了一个负号
利用这个特点，dft原本要算$N$个，现在只要算一半
根据$N$是$2$的整数幂，我们还可以继续分下去
$$\begin{array}{rl}& G(k)\\ & =\sum _{l=0}^{\frac{N}{4}-1}x(4l)W_{\frac{N}{2}}^{2lk}+\sum _{l=0}^{\frac{N}{4}-1}x(2(2l+1))W_{\frac{N}{2}}^{(2l+1)k}\\ & =\sum _{l=0}^{\frac{N}{4}-1}x(4l)W_{\frac{N}{4}}^{lk}+W_{\frac{N}{2}}^k\sum _{l=0}^{\frac{N}{4}-1}x(4l+2)W_{\frac{N}{4}}^{lk}\\ & =GG(k)+W_{\frac{N}{2}}^{k}GH(k)\end{array}$$
$$\begin{array}{rl}& H(k)\\ & =\sum _{l=0}^{\frac{N}{4}-1}x(4l+1)W_{\frac{N}{4}}^{lk}+W_{\frac{N}{2}}^k\sum _{l=0}^{\frac{N}{4}-1}x(4l+3)W_{\frac{N}{4}}^{lk}\\ & =HG(k)+W_{\frac{N}{2}}^{k}HH(k)\end{array}$$
所以，
周期为$N$的，可以划分为奇偶两组，每组再划分奇偶…
时间复杂度就是$O(N\log N)$
当然快速傅立叶逆变换非常类似，把$W_N$换成他的共轭$\overline{W_N}$就行了
参考代码

#include <cstdio>
#include <cstring>
#include <cmath>
using namespace std;
const double PI=acos(-1);
const int N=(1<<21)+5;//比2n大的2的整数幂
//复数
class Complex{
public:
    double real,image;
    explicit Complex(const double& real=0,const double& image=0):real(real),image(image){}
};
//复数加法
Complex operator+(const Complex& a,const Complex& b){
    return Complex(a.real+b.real,a.image+b.image);
}
//复数减法
Complex operator-(const Complex& a,const Complex& b){
    return Complex(a.real-b.real,a.image-b.image);
}
//复数乘法
Complex operator*(const Complex& a,const Complex& b){
    return Complex(a.real*b.real-a.image*b.image,b.real*a.image+a.real*b.image);
}
Complex temp[N];
/**
 * 快速傅立叶变换
 * @param a 变换数组
 * @param n N
 * @param inv 1为傅立叶变换，-1为傅立叶逆变换
 */
void fft(Complex a[],int n,int inv){
    if(n==1)return;
    int mid=n>>1; //N/2
    //划分奇偶
    for(int i=0;i<mid;++i){
        temp[i]=a[i<<1];
        temp[i+mid]=a[i<<1|1];
    }
    memcpy(a,temp,sizeof(Complex)*n);
 
    fft(a, mid, inv);
    fft(a+mid, mid, inv);
    //exp(2 PI i /N) || exp(-2 PI i /N)
    Complex w_n(cos(2*PI/n), inv * sin(2 * PI / n));
    Complex w(1,0);
    for(int i=0;i<mid;++i){
        Complex w_h=w*a[i+mid];
        //X(k)=G(k)+W_N^k H(k)
        temp[i]=a[i]+w_h;
        //X(k+N/2)=G(k)-W_N^k H(k)
        temp[i+mid]=a[i]-w_h;
 
        w=w*w_n;
    }
    memcpy(a,temp,sizeof(Complex)*n);
}
char s[N];
Complex a[N];
int ans[N];
int main(){
    scanf("%s",s);
    int len1=strlen(s);
    //我这个实部，你可以继续输入虚部
    for(int i=0;i<len1;++i){
        a[i].real=s[i]-'0';
    }
 
    int tot=len1;
    //变成2的n次幂
    int n=1;
    while(n<tot){
        n<<=1;
    }
    //快速傅立叶变换
    fft(a,n,1);
    //快速傅立叶逆变换
    fft(a,n,-1);
 
    for(int i=0;i<n;++i){
        //我这里只有实部，你可以再搞个处理虚部的
        ans[i]+=a[i].real/n+0.5;
    }
   
    return 0;
}

蝴蝶操作

上面的代码是递归的，而我们可以把他改成迭代的
以$N=4$为例
首先我们计算$G(0),G(1)$
$G(0)=x(0)+W_2^{0}x(2)$
$G(1)=x(0)-W_2^{0}x(2)$
对应到递归代码里，应该是$4$分成$2$个$2$个，然后$2$个分为$1$个$1$个（$x(0),x(2)$,直接return),然后在这里计算fft(a,2,1)

类似的计算H

接着计算$X$
$X(0)=G(0)+W_4^{0}H(0)$
$X(1)=G(1)+W_4^{1}H(1)$
$X(2)=G(2)-W_4^{0}H(2)$
$X(3)=G(3)-W_4^{1}H(3)$

可以看出来，
由$x$计算$G$和$H$的时候
  旋转因子的$N$为$2$
  可以分为$2$组，组里下标相差$\frac{N}{4}$的加在了一起
而$G$和$H$计算$X$的时候
  旋转因子的$N$为$4$
  可以分为$1$组，下标相差$\frac{N}{2}$的加在了一起
旋转因子每次只会循环到$N/2$，因为另外$N/2$都是利用对称性的
最左边的$x$是根据$1$次划分奇偶排的

我们接着看$N=8$的
可以看出来也差不多
$x$计算$GG,GH,HG,HH$
  旋转因子的$N$为$2$
  可以分为$4$组，组里下标相差的加在了一起
$GG,GH,HG,HH$计算$G，H$
  旋转因子的$N$为$4$
  可以分为2组，组里下标相差的加在了一起
$G,H$计算$X$
  旋转因子的$N$为$8$
  可以分为$1$组，组里下标相差的加在了一起
旋转因子每次只会循环到$N/2$，因为另外$N/2$都是利用对称性的
最左边的$x$是根据$2$次划分奇偶排的
0 1 2 3 4 5 6 7
0 2 4 6 1 3 5 7
0 4 2 6 1 5 3 7

先解决这个奇偶划分
  0  1  2  3    4   5   6  7
000 001 010 011 100 101 110 111
反过来
  0  1  2  3    4   5   6  7
000 100 010 110 001 101 011 111
  0  4  2  6    1   5   3  7
于是，我们可以写出

rev[i]=rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));

接着就是观察规律了
于是不难写出代码

#include <cstdio>
#include <cstring>
#include <cmath>
#include<iostream>
using namespace std;
const double PI=acos(-1);
const int N=(1<<21)+5;//比2n大的2的整数幂
//复数
class Complex{
public:
    double real,image;
    explicit Complex(const double& real=0,const double& image=0):real(real),image(image){}
};
//复数加法
Complex operator+(const Complex& a,const Complex& b){
    return Complex(a.real+b.real,a.image+b.image);
}
//复数减法
Complex operator-(const Complex& a,const Complex& b){
    return Complex(a.real-b.real,a.image-b.image);
}
//复数乘法
Complex operator*(const Complex& a,const Complex& b){
    return Complex(a.real*b.real-a.image*b.image,b.real*a.image+a.real*b.image);
}
 
int rev[N];
/**
 * 快速傅立叶变换
 * @param a 变换数组
 * @param n N
 * @param inv 1为傅立叶变换，-1为傅立叶逆变换
 */
void fft(Complex a[],int n,int inv){
    for(int i=0;i<n;++i){
        if(i<rev[i]){
            swap(a[i],a[rev[i]]);
        }
    }
    //i=每层的N  /2   (也就是每组要循环多少下）
    for(int i=1;i<n;i<<=1){
        //exp(-2 PI i /N) || exp(2 PI i /N)
        Complex w_n(cos(PI/i),inv*sin(PI/i));
        //j是每组开始下标
        for(int j=0;j<n;j+=(i<<1)){
            Complex w(1,0);
            //k是组里的第几个
            for(int k=0;k<i;++k,w=w*w_n){
                Complex g=a[j+k],w_h=w*a[j+k+i];
                a[j+k]=g+w_h;
                a[j+k+i]=g-w_h;
            }
        }
    }
 
    if(-1 == inv){
        for(int i=0;i<n;++i){
            a[i].real=(int)(a[i].real/n+0.5);
            a[i].image=(int)(a[i].image/n+0.5);
        }
    }
}
char s[N];
Complex a[N];
int main(){
    scanf("%s",s);
    int len1=strlen(s);
    //我这个实部，你可以继续输入虚部
    for(int i=0;i<len1;++i){
        a[i].real=s[i]-'0';
    }
 
    int tot=len1;
    int n=1,len=0;
    while(n<tot){
        n<<=1;
        ++len;
    }
    for(int i=0;i<n;++i){
        rev[i]=rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
    }
    fft(a,n,1);
    //快速傅立叶逆变换
    fft(a,n,-1);
    
    return 0;
}

大数乘法

其实我大一的时候，写大数乘法，不会写，然后就去搜，搜到了fft，给我整不会了
众所周知
$$\mathfrak{F}\left[f\ast g\right]=\mathfrak{F}\left[f\right]\mathfrak{F}\left[g\right]$$
然后大数乘法的每个系数其实就是卷积，于是
我们可以读入两个大数，然后逆序存储
对他们分别傅立叶变换，然后对应位相乘，最后逆变换回来，就是最终的答案了
注意一下，$n$位数$\ast m$位数，不超过$n+m+1$位数，然后你要把$n+m+1$补到2的幂
以洛谷P1919为例

#include <cstdio>
#include <cstring>
#include <cmath>
#include<iostream>
using namespace std;
const double PI = acos(-1);
const int N = (1 << 21) + 5;//比2n大的2的整数幂
//复数
class Complex {
public:
    double real, image;
    explicit Complex(const double& real = 0, const double& image = 0) :real(real), image(image) {}
};
//复数加法
Complex operator+(const Complex& a, const Complex& b) {
    return Complex(a.real + b.real, a.image + b.image);
}
//复数减法
Complex operator-(const Complex& a, const Complex& b) {
    return Complex(a.real - b.real, a.image - b.image);
}
//复数乘法
Complex operator*(const Complex& a, const Complex& b) {
    return Complex(a.real * b.real - a.image * b.image, b.real * a.image + a.real * b.image);
}
 
int rev[N];
/**
 * 快速傅立叶变换
 * @param a 变换数组
 * @param n N
 * @param inv 1为傅立叶变换，-1为傅立叶逆变换
 */
void fft(Complex a[], int n, int inv) {
    for (int i = 0; i < n; ++i) {
        if (i < rev[i]) {
            swap(a[i], a[rev[i]]);
        }
    }
    //i=每层的N  /2   (也就是每组要循环多少下）
    for (int i = 1; i < n; i <<= 1) {
        //exp(-2 PI i /N) || exp(2 PI i /N)
        Complex w_n(cos(PI / i), inv * sin(PI / i));
        //j是每组开始下标
        for (int j = 0; j < n; j += (i << 1)) {
            Complex w(1, 0);
            //k是组里的第几个
            for (int k = 0; k < i; ++k, w = w * w_n) {
                Complex g = a[j + k], w_h = w * a[j + k + i];
                a[j + k] = g + w_h;
                a[j + k + i] = g - w_h;
            }
        }
    }
 
    if (-1 == inv) {
        for (int i = 0; i < n; ++i) {
            a[i].real = (int)(a[i].real / n + 0.5);
            //a[i].image=(int)(a[i].image/n+0.5);
        }
    }
}
char s[N];
Complex a[N], b[N];
int ans[N];
 
int main() {
    scanf("%s", s);
    int len1 = 0;
    for (int i = strlen(s) - 1; i >= 0; --i) {
        a[len1++].real = s[i] - '0';
    }
    scanf("%s", s);
    int len2 = 0;
    for (int i = strlen(s) - 1; i >= 0; --i) {
        b[len2++].real = s[i] - '0';
    }
    int tot = len1 + len2 + 1;
    int n = 1, len = 0;
    while (n < tot) {
        n <<= 1;
        ++len;
    }
    for (int i = 0; i < n; ++i) {
        rev[i] = rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len - 1));
    }
    fft(a, n, 1);
    fft(b, n, 1);
    for (int i = 0; i < n; ++i) {
        a[i] = a[i] * b[i];
    }
    fft(a, n, -1);
 
    for (int i = 0; i < n; ++i) {
        ans[i] += a[i].real;
        ans[i + 1] += ans[i] / 10;
        ans[i] %= 10;
    }
    ++n;
    while (n > 0 && 0 == ans[n])--n;
    while (n >= 0) {
        printf("%d", ans[n]);
        --n;
    }
    printf("\n");
    return 0;
}