Fork me on GitHub

FFT残疾人手册

Posted on | Edited on | In

FFT残疾人手册

0.简介

FFT/NTT可以把本来复杂度是$\mathcal O(n^2)$的多项式乘法优化成$\mathcal O(n\log n)$的复杂度

在比较高等级的比赛中会用吧

在电脑上存一个板子也挺方便的

1.点值表示法

一个$n$次多项式

的点值表示法为

其实就是把$x_0,x_1,…,x_{n-1}$代入$A(x)$中得到的$y_0,y_1,…,y_{n-1}$

一个大小为$n$点值表示法可以唯一确定一个$n$次多项式

因为点值表示法实际上就是n个方程

这样是可以唯一确定系数集合$a_0,a_1,…,a_{n-1}$的

2.算法的主要思想

我们把多项式$A(x),B(x)$在$x_0,x_1,…,x_{2n-2}$处的点值分别求出来，得到

把这些点值乘起来变成

就可以还原成多项式$A(x)*B(x)$

3.前置知识：单位根

给没有复数基础知识的OIer普及一下：$i=\sqrt{-1}$，是一个常数

然后形如$cos(\theta)+i\cdot sin(\theta)$的复数有一个很好的性质

$(cos(\theta)+i\cdot sin(\theta))^k=cos(k\theta)+i\cdot sin(k\theta)$

这个东西可以用数学归纳法暴力展开证明

n次单位根是形如$cos(\frac{2\pi k}{n})+i\cdot sin(\frac{2\pi k}{n})$的复数

我们记$\omega_{n}=cos(\frac{2\pi}{n})+i\cdot sin(\frac{2\pi}{n})$

然后有$\omega_n^k=cos(\frac{2\pi k}{n})+i\cdot sin(\frac{2\pi k}{n})$

所以说n次单位根其实就是把单位元分成了n等份

根据这个性质很容易得到

也就是说存在：

根据三角函数的周期性可以得到一个式子：

结合$(7)(8)$可以得到：

还有一个很好的性质：

第$(9)$和$(10)$个式子都是我们要用到重要性质

4.DFT过程及其优化

DFT就是所谓的离散傅里叶变换（Discrete Fourier Transform, DFT）

因为朴素的多项式转化成点值的复杂度是$n^2$的，很慢

所以考虑怎么把这个东西优化下来

首先DFT的优化用到了分治的思想

根据第$(9)$个式子，我们来猜想一些性质

多项式：

把A的所有项按照x指数的奇偶性来分类变成

令：

所以有：

那么$A(x)$就变成了$F(x)$在$x^2$处的点值和$G(x)$在$x^2$处的点值和x的积之和

在取x集合的时候令：

就可以得到：

这样的话我们只需要计算出$F(x)$和$G(x)$的值就可以算出多项式$A(x)$的点值了

而因为每次我们递归的时候都会缩小一半的计算范围，只不过需要预先把n补到2的次幂

所以这样的时间复杂度是：

5.IDFT过程及其优化

我们做完了DFT，现在需要用点值还原出多项式

也就是说要从这个有n个等式的线性方程组：

解出$a_0,a_1,…a_{n-1}$的值。

我们把上面的线性方程组写成矩阵形式，发现是：

也就是说我们要求出矩阵

的逆矩阵$P’$。

在这里直接给出逆矩阵$P’$并证明一定成立

那么对于矩阵$Q=PP’$，有$Q_{i,j}=\sum_{k=0}^{n-1}P_{i,k}\cdot P’_{k,j}$

即$Q_{i,j}=\sum_{k=0}^{n-1}\frac{1}{n}\omega_n^{-ik}\cdot \omega_n^{kj}=\frac{1}{n}\sum_{k=0}^{n-1}(\omega_n^{j-i})^k $

  1. $i\not=j$：$Q_{i,j}=\frac{1}{n}\sum_{k=0}^{n-1}(\omega_n^{i-j})^k=\frac{1}{n}\cdot \frac{(\omega_n^{i-j})^n-1}{\omega_n^{i-j}-1}$

    因为公式$(8)$，可以得到$(\omega_n^{i-j})^n=1$，所以$Q_{i,j}=0$

  2. $i=j$：$Q_{i,j}=1$

所以$P​$和$P’​$互为逆矩阵

也就是说现在要做的就是求出$P’$和点值矩阵的乘积，即：

再来观察一下原来的DFT过程中我们做的事，就是求出了

那这个时候我们只需要把$A(\omega_n^0) ,A(\omega_n^1) ,\dots,A(\omega_n^{n-1})$当作点值，并且把DFT过程中的$\omega_n^{k}$替换成$\omega_n^{-k}$就可以了

6.FFT的递归实现

dd只要理解了最初的迭代过程应该就问题不大了

就是简单的模拟出分治过程

omiga数组最好是可以预处理出来，但是这样的实现方式依旧有很大的常数

1
2
3
4
5
6
7
8
9
10
11
12
void fft(int n, Complex *p, int pos, int step, Complex *omg) {
  if(n == 1) return;
  int m = n >> 1;
  fft(m, p, pos, step << 1, omg);
  fft(m, p, pos + step, step << 1, omg);
  for(int k = 0; k < m; k++) {
    int cur = 2 * step * k;
    tmp[k] = p[cur + pos] + omg[k * step] * p[cur + pos + step];
    tmp[k + m] = p[cur + pos] - omg[k * step] * p[cur + pos + step];
  }
  for(int i = 0; i < n; i++) p[i * step + pos] = tmp[i];
}

7.FFT的迭代实现

我们考虑每一次的分组递归

第一次显然是按照最后一个二进制位的奇偶性来分组的

进入下一层迭代后就发现变成了一个$\frac{n}{2}$次多项式，可以看做所有的指数都除以了2，参考$(14)$

所以最后一位就没有意义了，下一次排序就是按照倒数第二位的奇偶性来进行分组的

发现本质就是把二进制反转进行交换，很奇妙的

可以这样预处理：

1
2
3
for (int i = 0; i < n; i++) {
  rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (lim - 1));
}

所以就可以用迭代来进行运算啦

1
2
3
4
5
6
7
8
9
10
11
12
13
void transform(Complex f[N], Complex w[N]) {
  for (int i = 0; i < p; i++) if (i < rev[i]) swap(f[i], f[rev[i]]);
  for (int len = 1; len < p; len <<= 1) {
    int step = (len << 1);
    for (Complex *cur = f; cur != f + p; cur += step) {
      for (int i = 0; i < len; i++) {
        Complex t = w[p / step * i] * cur[i + len];
        cur[i + len] = cur[i] - t;
        cur[i] += t; 
      }
    }
  }
}

8.板子

总的板子

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
#include<bits/stdc++.h>

using namespace std;

const double PI = acos(-1.0);
const int N = 262144;

struct Complex {
	double x, y;
	
  Complex() {}
	
  Complex(double x, double y): x(x), y(y) {}
	
  Complex operator + (const Complex & a) const { 
    return Complex(x + a.x, y + a.y);
  }
	
  void operator += (const Complex & a) {
    *this = *this + a;
  }
	
  Complex operator - (const Complex & a) const {
    return Complex(x - a.x, y - a.y);
  }
  
  void operator -= (const Complex & a) {
    *this = *this - a;
  }
  
  Complex operator * (const Complex & a) const {
    return Complex(x * a.x - y * a.y, x * a.y + y * a.x);
  }
  
  void operator *= (const Complex & a) {
    *this = *this * a;
  }
  
  Complex operator / (const Complex & b) const {
    Complex c = b.conj();
    return (*this) * c / ((b * c).x);
  }
  
  void operator /= (const Complex & a) {
    *this = *this / a;
  }
  
	Complex operator * (const double & a) const { 
    return Complex(x * a, y * a);
  }
  
  void operator *= (const double & a) {
    *this = *this * a;
  }
  
	Complex operator / (const double & a) const {
    return Complex(x / a, y / a);
  }
  
  void operator /= (const double & a) {
    *this = *this / a;
  }
  
	Complex conj() const {
    return Complex(x, -y);
  }
} ;

struct fft {
  Complex a[N], b[N], c[N], omg[N], inv[N];
  int n, m, p, lim, rev[N];
  
  void input() {
   scanf("%d %d", &n, &m);
    ++n, ++m;
    for (int i = 0; i < n; i++) {
      int v; scanf("%d", &v);
      a[i] = Complex(v, 0);
    }
    for (int i = 0; i < m; i++) {
      int v; scanf("%d", &v);
      b[i] = Complex(v, 0);
    }
  }
  
  void output() {
    for (int i = 0; i < n + m - 1; i++) {
      printf("%d ", (int)(c[i].x + 0.5));
    }
  }
  
  void init() {
    for (p = 1; p < n + m - 1; p <<= 1);
    for (; (1 << lim) < p; ++lim);
    for (int i = 0; i < p; i++) {
      omg[i] = Complex(cos(2.0 * PI * i / p), sin(2.0 * PI * i / p));
      inv[i] = omg[i].conj();
      rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (lim - 1));
    }
  }
  
  void transform(Complex f[N], Complex w[N]) {
    for (int i = 0; i < p; i++) if (i < rev[i]) swap(f[i], f[rev[i]]);
    for (int len = 1; len < p; len <<= 1) {
      int step = (len << 1);
      for (Complex *cur = f; cur != f + p; cur += step) {
        for (int i = 0; i < len; i++) {
          Complex t = w[p / step * i] * cur[i + len];
          cur[i + len] = cur[i] - t;
          cur[i] += t; 
        }
      }
    }
  }
  
  void solve() {
    input();
    init();
    transform(a, omg);
    transform(b, omg);
    for (int i = 0; i < p; i++) c[i] = a[i] * b[i];
    transform(c, inv);
    for (int i = 0; i < p; i++) c[i] = Complex(c[i].x / p, 0.0);
    output();
  }
} fft;

int main() {
  fft.solve();
  return 0;
}