【数论系列】快速傅里叶变换 FFT算法

阿阿阿安

已于 2024-06-07 21:02:54 修改

阅读量5.3k

点赞数 6

分类专栏：【算法竞赛】 # 数学几何文章标签：算法

于 2019-07-25 21:51:14 首次发布

本文链接：https://blog.csdn.net/qq_40772692/article/details/97306446

版权

【算法竞赛】同时被 2 个专栏收录

186 篇文章 23 订阅

订阅专栏

数学几何

52 篇文章 2 订阅

订阅专栏

3. FFT快速傅里叶变换与快速傅里叶逆变换

2. 计算DFT(A(x))和DFT(B(x))

一. FFT基础

1. 基本概念

定义可以自行百度。通俗点来说，FFT就是利用某些奇偶特点，进行DFT（离散傅里叶变换）和IDFT（离散傅里叶逆变换）的快速求解算法。

2. 应用场景

（1）该算法在信号学中有很大用处；

（2）在信息学竞赛中：加速多项式乘法，高精度大数运算等；

二. FFT的发展历程

已知多项式:

$A(x) = a0 + a1*x + a2*x^{2} + ....+a(n-1)*x^{n-1}$

$B(x) = b0 + b1*x + b2*x^{2} + ....+b(n-1)*x^{n-1}$

现在要计算两个n项的n-1次多项式相乘，其结果肯定是2n-1项的2n-2次多项式，求出相乘多项式幂次由低到高的系数。

1. 朴素算法阶段（暴力）

普通的多项式乘法是O(n^2)的，我们要枚举A(X)中的每一项，分别与B(x)中的每一项相乘，来得到一个新的多项式C(x)。

2. 傅里叶的DFT和IDFT

2.1 多项式的表示方法

（1）系数表示法：用幂次不定元及其系数表示的多项式，比如A（X）。

（2）点值表示法：任何一个n次多项式，都能由n+1个不同的点唯一确定，用这n+1个点表示的多项式称为多项式的点值表示。

2.2 DFT和IDFT的思路

已知A(X)和B(X)的系数表示，求C(X) = A(X)*B(X)；

我们将A(X)与B(X)先转化为点值表示法，即：A(X) = { (p0,A(p0))， (p1,A(p1))，(p2,A(p2))，.....} 、B(X) = { (p0,B(p0))， (p1,B(p1))，(p2,B(p2))，..... }；

则在点值表示法下：计算点值C(pi) = A(pi)*B(pi) 是 O(n) 的，这个过程叫做DFT。然后我们再通过IDFT将C(pi)的点值表示法转化为系数表示法，这个过程插值。但是遗憾的是：系数表示转化为点值表示是O(n^2) ，插值过程也是O(n^2) 。在计算机中并不能提高效率。

3. FFT快速傅里叶变换与快速傅里叶逆变换

（1）目的：我们想要利用DFT中点值表示法计算的O(n)复杂度，又想缩短转化过程中的时间复杂度。于是就有了FFT算法，其本质是利用分治的思想加快DFT和IDFT的计算过程，通过各种优化使得整个求解过程降为 O(nlogn)。

（2）过程：见下文。

三. 前置知识

1. 复数基础

1.1 基本性质

（1）定义：设a,b为实数， $i^{2} = -1$ ，形如 a + bi 的数叫复数，其中 i 被称为虚数单位。在复平面中，x轴代表实数，y轴（除原点外的点）代表虚数，从原点(0,0)到(a,b)的向量表示复数 a + bi ；

（2）模长：从原点(0,0)到点(a,b)的距离，即 $\sqrt{a^{2} + b^{2}}$ ；

（3）幅角：假设以逆时针为正方向，从x轴正半轴到已知向量的转角的有向角叫做幅角；

1.2 运算规则

（1）几何定义：复数相乘，模长相乘，幅角相加（所以说模长为1的复数相乘永远模长为1）；

（2）代数定义：(a+bi)∗(c+di) = (ac−bd)+(bc+ad)i；

1.3 单位根

（1）基本概念

已知： $z^{n} = 1(n = 1,2,3...)$ ；则满足这个方程的复数解z ，称为n次单位根，记为 $w_{n}$ ；这样的的 n次根有n个：它们位于复平面的单位圆上，构成正n边形的顶点，其中一个顶点是1。

根据复数乘法的运算法则，其余n个复数为 $w_{n}^{1},w_{n}^{2},w_{n}^{3}.....w_{n}^{n}(w_{n}^{0} = w_{n}^{n} = 1)$ 且 $w_{n}^{k} = (w_{n}^{1})^k$ ，那么如何计算它们的值呢？这个问题可以由欧拉公式解决：

如八次单位根如下：

（2）基本性质（FFT就是利用这个性质简化运算）

性质1： $w_{2n}^{2k} = w_{n}^{k}$ ；
性质2： $w_{n}^{k + \frac{n}{2}} = -w_{n}^{k}$ ；

四. 快速傅里叶变换

1. 前置条件

（1）已知：一个n-1次多项式，可以被n个不同的点值唯一确定。

（2）所以在求点值过程中，我们必须选择n个x带入方程，来得到n个点值为该方程的点值表示。

（3）那么这n个自变量x怎么选择呢？前面我们讲了这么多，当然是代入n次单位根的n个值 $w_{n}^{0},w_{n}^{1},w_{n}^{2}.....w_{n}^{n-1}$

（4）接下来是利用单位根的性质来简化运算，加速求值过程。

2. 推导过程

设所要求值多项式为： $A(x) = a0 + a1*x + a2*x^{2} + ....+a(n-1)*x^{n-1}$ （n必须均为2的幂次，原因在下文），则：

$A(x) = (a0 + a2*x^{2} + a4*x^{4} + ..... + a_{n-2}*x^{n-2})+(a1*x + a3*x^{3} + .... + a_{n-1}*x^{n-1})$

$A(x) = (a0 + a2*x^{2} + a4*x^{4} + ..... + a_{n-2}*x^{n-2})+x*(a1 + a3*x^{2} + .... + a_{n-1}*x^{n-2})$

令 $A_{1}(x) = a0 + a2*x+ a4*x^{2} + ..... + a_{n-2}*x^{\frac{n}{2}-1}$

令 $A_{2}(x) = a1 + a3*x+ a5*x^{2} + ..... + a_{n-1}*x^{\frac{n}{2}-1}$

则 $A(x) = A_{1}(x^{2})+ x*A_{2}(x^{2})$ ，那么：

（1）假设我们代入 $w_{n}^{k}$ 时 (k < n/2)：性质1

$A(w_{n}^{k}) = A_{1}(w_{n}^{2k})+ w_{n}^{k}*A_{2}(w_{n}^{2k}) = A_{1}(w_{\frac{n}{2}}^{k})+ w_{n}^{k}*A_{2}(w_{\frac{n}{2}}^{k})$ ；

（2）同理当代入另一半值时代入 $w_{n}^{k + \frac{n}{2}}$ (k<n/2)：性质1 + 性质2 + $w_{n}^{n} = 1$ $A(w_{n}^{k+\frac{n}{2}}) = A_{1}(w_{n}^{2k+n})+ w_{n}^{k+\frac{n}{2}}*A_{2}(w_{n}^{2k+n}) = A_{1}(w_{\frac{n}{2}}^{k} * w_{n}^{n}) - w_{n}^{k}*A_{2}(w_{\frac{n}{2}}^{k} * w_{n}^{n}) = A_{1}(w_{\frac{n}{2}}^{k}) - w_{n}^{k}*A_{2}(w_{\frac{n}{2}}^{k})$

综上所述：我们只需要n/2个值，就能求出A(x)的n个点值表示，可以看出来，我们将问题的规模缩小了一半！只要我们不断往下分治，最后到常数项返回，复杂度会缩短为O(nlogn)！

五. 快速傅里叶逆变换

1. 推导过程

参考文章：从多项式乘法到快速傅里叶变换 | Miskcoo’s Space

（1）有了最终多项式的点值表示法以后，我们要把点值表示法转化为系数表示法：

（2）假设方程为 W*a = A，则方程左侧同时乘以 $W^{-1}$ （W逆）得： $a = W^{-1} * A$ ，则经过计算得：

即快速IDFT的过程是以所求点值表示的A(x)作为某方程T(x)系数，代入n个 $x = w_{n}^{-k}$ 所求得的T(x)的点值表示中，T(xi) / n 即为所转化系数。也就是用 $x = w_{n}^{-k}$ 再进行一次快速傅里叶变换的过程。

注意：对于模长为1的复数，其倒数等同于其共轭复数。即 $w_{n}^{-k} = a-bi$ ( $w_{n}^{k} = a+bi$ )；

六. FFT实现过程

1. 补项操作

已知n+1项n次多项式A(X)和m+1项m次多项式B(X),求C(X) = A(X)*B(X)。

则C（X）最多为n+m+1项n+m次多项式，我们需要将A(X)和B(X)均补齐到>=n+m+1的最小2幂次项limit，缺的幂次前面系数补0（为了分治过程中满足每次都能均分为奇偶分段），在输出时不输出多余项即可。

2. 计算DFT(A(x))和DFT(B(x))

需要代入limit次单位根的limit个值,即 $w_{limit}^{k} = \cos (k*\frac{2\pi }{limit}) + i*\sin(k*\frac{2\pi }{limit})$ ，而STL提供了complex复数模板，complex.real( )为实数部分，complex.imag()为虚数部分（也可手动高精度实现），分离奇偶，递归求出 $A(w_{\frac{n}{2}}^{k})$ 然后求出 $A(w_{n}^{k})$ 和 $A(w_{n}^{k + n/2})$ 。

3. 求C(x)的点值表示

即C(X) = A(X)*B(X)。

4. 计算IDFT(C(X))

需要代入limit次单位根倒数的limit个值,即 $w_{limit}^{-k} = \cos (k*\frac{2\pi }{limit}) - i*\sin(k*\frac{2\pi }{limit})$

，注意结果不要忘记 /n。

七. 代码实现及优化

洛谷P3803 【模板】多项式乘法（FFT）

1. 基础递归版算法

#include<iostream>
#include<bits/stdc++.h>
#define PI acos(-1)
using namespace std;
const int maxn = 1000000 + 7;
int n,m;
complex<double> a[maxn*3],b[maxn*3];
void FFT(int len,int type,complex<double> *c){
   if(len==1)return;
   complex<double> codd[len>>1];
   complex<double> ceven[len>>1];
   for(int i = 0;i<len;i+=2){//奇偶分离A(X) = A1(X) + x*A2(X)
       ceven[i>>1] = c[i];//A1(X)
       codd[i>>1] = c[i+1];//A2(X)
   }
   FFT(len>>1,type,ceven);//求子问题
   FFT(len>>1,type,codd);
   for(int i = 0;i<(len>>1);i++){//套公式
      complex<double> wn(cos(2.0*PI*i/len),type*sin(2.0*PI*i/len));
      c[i] = ceven[i] + wn*codd[i];
      c[i+(len>>1)] = ceven[i] - wn*codd[i];
   }

}

int main(){
  int x;
  int maxLen = 1;
  scanf("%d%d",&n,&m);
  while(maxLen<n+m+1)maxLen<<=1;//补项到最终结果最小的2的幂次项
  for(int i = 0;i<=n;i++){
      scanf("%d",&x);
      a[i].real(x);//n+1个A(X)系数a0......an
  }
  for(int i = 0;i<=m;i++){
      scanf("%d",&x);
      b[i].real(x);//m+1个B(X)系数b0....bm
  }

  FFT(maxLen,1,a);//DFT求点值（nlogn）
  FFT(maxLen,1,b);

  for(int i = 0;i<maxLen;i++){//求C(X)点值
     a[i]*=b[i];
  }
  FFT(maxLen,-1,a);//IDFT
  for(int i = 0;i<=n+m;i++){//只输出到n+m+1项
    if(i!=0)printf(" ");
    printf("%d",(int)(a[i].real()/maxLen+0.5));
  }
  printf("\n");
  return 0;
}

2. 迭代求解+蝴蝶操作优化

首先找一下每个数字递归到最后所在的位置关系：

我们发现：每一个位置处最后的位置是该位置最初二进制的反转，所以我们怎么实现呢？在原序列中：

（1）若位置i为偶数：位置i的反转位置 = i/2的二进制左移一位

（2）若位置i为奇数：位置i的反转位置 = i/2的二进制左移一位后在最高位补1

所以我们用数组保存下每个原序列位置最后的位置：pos[i] = (pos[i>>1]>>1)|((i&1)<<(l-1)) 。交换完毕以后，这就相当于递归到的最后一层，我们需要不断向上合并得到最终答案。

注意：使用complex<double> Wn (cos(2.0*PI/L),type*sin(2.0*PI/L)); + complex<double>w(1,0); + w=w*Wn的方式来求 $W_{n}^{k}$ 可以使时间缩短一半，预处理的话会更快。

#include <iostream>
#include<bits/stdc++.h>
using namespace std;
const int maxn = 1000000 + 7;
#define PI acos(-1)
int n,m;
complex<double> a[maxn*3],b[maxn*3];
int pos[maxn*3];
void FFT(complex<double>*A,int len,int type){
   for(int i = 0;i<len;i++){//把每个数放到最后的位置
       if(i<pos[i])swap(A[i],A[pos[i]]);//保证每对只交换一次
   }
   for(int L = 2;L<=len;L<<=1){//循环合并的区间长度
     int HLen = L/2;//区间的一半
     complex<double> Wn (cos(2.0*PI/L),type*sin(2.0*PI/L));
     for(int R = 0;R<len;R+=L){//每个小区间的起点
        complex<double>w(1,0);
        for(int k = 0;k<HLen;k++,w=w*Wn){//求该区间下的值
            complex<double> Buf = A[R+k];//蝴蝶操作，去掉odd和even数组，使变化原地进行
            A[R+k] =  A[R+k] + w*A[R+k+HLen];
            A[R+k+HLen] = Buf - w*A[R+k+HLen];
        }
     }
   }
}
int main()
{
    int x;pos[0] = 0;
    int maxLen = 1,l = 0;
    scanf("%d%d",&n,&m);
    while(maxLen<n+m+1){maxLen<<=1;l++;}
    for(int i = 0;i<=n;i++){
        scanf("%d",&x);a[i].real(x);
    }
    for(int i = 0;i<=m;i++){
        scanf("%d",&x);b[i].real(x);
    }
    for(int i = 0;i<maxLen;i++){//求最后交换的位置（奇数特殊处理）
        pos[i] = (pos[i>>1]>>1)|((i&1)<<(l-1));
    }
    FFT(a,maxLen,1);
    FFT(b,maxLen,1);
    for(int i = 0;i<maxLen;i++)a[i]*=b[i];
    FFT(a,maxLen,-1);
    for(int i = 0;i<n+m+1;i++){
        if(i!=0)printf(" ");
        printf("%d",(int)(a[i].real()/maxLen+0.5));
    }
    printf("\n");
    return 0;
}