FFT学习笔记
快速傅里叶变换
多项式
定义
不严谨定义：形如 \(f(x) = \sum \limits _{i=0}^{n} a_ix^i\) 的式子为多项式。
定义（from OIWiki）：对于求和式 \(\sum a_nx^n\)，如果是有限项相加，称为多项式，记作 \(f(x)=\sum \limits_{n=0}^m a_nx^n\)。
次数：对于多项式 \(F(x) = \sum \limits _{i=0}^{n} a_ix^i\)，该多项式的次数为 \(n\)。表示为 \(\text{degree} (F) = n\)
次数界：任何一个大于多项式次数的整数。
基本运算
系数表示法：
加法
时间复杂度 \(O(n)\)。
乘法
记为 \(\mathbf{c} = \mathbf{a} \otimes \mathbf{b}\)
时间复杂度 \(O(n^2)\)。
求值
秦九韶算法：
时间复杂度 \(O(n)\)。
点值表示法：
加法
时间复杂度 \(O(n)\)。
乘法
时间复杂度 \(O(n)\)。
求值
拉格朗日插值法：
时间复杂度 \(O(n^2)\)
代码实现：
double res = 0;
for(int i = 1; i <= n; i ++ )
{
double s1 = 1, s2 = 1;
for(int j = 1; j <= n; j ++ )
if(j != i)
{
s1 = s1 * (k - x[j]);
s2 = s2 * (x[i] - x[j]);
}
res = res + y[i] * s1 / s2;
}
复数
基本运算
高中数学基础，不多介绍
欧拉定理
欧拉定理：\(e^{i\theta} = \cos \theta + i\sin \theta\)
根据这个定理，我们可以把复数转化成在复平面上的一堆点来进行一些公式的推导。
算术基本定理
\(n\) 次代数方程 \(f(x)=a_nx^n+a_{n-1}x^{n-1}+\cdots+a_1x+a_0 = 0(a_n \ne 0 且 n > 0)\) 在复数域上恰好有 \(n\) 个根。
根据这个定理，我们可以找出一个特殊的式子：\(x^n=1\) 的根，在复数域中恰好有 \(n\) 个，设 \(\omega _n = e^{\frac{2\pi}{n}i}\)，则 \(x^n = 1\) 的解集可以记作 \(\{ \omega _n^k|k=0,1,2,\cdots,n-1 \}\)。
单位根
单位根具有一些性质，在 FFT 的时候会有很大用处。
折半引理
这个结论在复平面中很好证明，当然也可以用公式来证明：
消去引理
证明：
求和引理
证明：
FFT
在计算多项式乘法时，系数表示法的复杂度为 \(O(n^2)\)，而点值表示法的复杂度仅为 \(O(n)\)，我们可以尝试将系数表示法转化为点值表示法来进行运算。
我们可以对多项式进行如下变换，将奇数项和偶数项分开：
（前提：\(n\) 为奇数且 \(n+1\) 为 \(2\) 的幂）
将奇数项的系数构成的多项式记为 \(A^{[1]}\)，偶数项记为 \(A^{[0]}\)，就可以得到下面的式子。
这时候我们想到，可以对两边的多项式分治处理，但是实数中没有可以把这个平方消掉的数，但是我们可以从复数中找到。
将 \(\omega _{2n}^k\) 带入 \(A\) 中，就可以得到：
我们就可以对 \(A^{[0]}\) 和 \(A^{[1]}\) 进行递归处理。
同时，在复平面中，\(\omega _{2n}^{n+k}\) 与 \(\omega _{2n}^{k}\) 在多项式 \(A^{[0]}(x^2)\) 和 \(A^{[1]}(x^2)\) 的值相同，因此我们也可以算出该点的值。
这就是 FFT 的实现原理。做完 FFT 之后，我们就得到了一个多项式在复数域中 \(2n\) 个点的点值，这时再将两个多项式乘起来，就可以得到乘积的点值表示。
但是如果我们需要的是系数表示，我们则需要做一遍快速傅里叶逆变换，将点值转化为系数。
我们可以列出刚才求点值时的矩阵表示：
我们现在只需要求出逆矩阵即可。由于这个矩阵元素非常特殊（没记错应该叫范德蒙德矩阵），它的逆矩阵即为每一项的倒数再除以变换的长度。
为了取倒数，我们可以做变换：
我们可以看到与 FFT 时仅差了一个负号，因此我们可以将 FFT 和 IDFT 结合为一题，传参时传入一个 \(opt\)，为 \(1\) 时执行 FFT，为 \(-1\) 时执行 IDFT。
代码实现：
#include <bits/stdc++.h>
using namespace std;
const double PI = acos(-1);
const int N = 3e6 + 10;
struct Complex
{
double a, b;
Complex() {a = 0, b = 0;}
Complex(double real, double imag): a(real), b(imag) {}
Complex operator + (const Complex& x) const
{
return Complex(a + x.a, b + x.b);
}
Complex operator - (const Complex& x) const
{
return Complex(a - x.a, b - x.b);
}
Complex operator * (const Complex& x) const
{
return Complex(a * x.a - b * x.b, a * x.b + b * x.a);
}
}F[N], G[N];
void FFT(Complex a[], int lim, int opt)
{
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(a1, lim >> 1, opt);
FFT(a2, lim >> 1, opt);
Complex wn = Complex(cos(2.0 * PI / lim), opt * sin(2.0 * PI / lim));
Complex w = Complex(1, 0);
for(int k = 0; k < (lim >> 1); k ++ )
{
a[k] = a1[k] + w * a2[k];
a[k + lim / 2] = a1[k] - w * a2[k];
w = w * wn;
}
}
int n, m;
int main()
{
n = read(), m = read();
for(int i = 0; i <= n; i ++ ) scanf("%lf", &F[i].a);
for(int i = 0; i <= m; i ++ ) scanf("%lf", &G[i].a);
int lim = 1;
while(lim <= n + m) lim <<= 1;
FFT(F, lim, 1), FFT(G, lim, 1);
for(int i = 0; i <= lim; i ++ ) F[i] = F[i] * G[i];
FFT(F, lim, -1);
for(int i = 0; i <= n + m; i ++ ) cout << (int)(F[i].a / lim + 0.5) << ' ';
return 0;
}
优化
BitReverse
上面我们使用了递归实现 FFT，我们可以思考一下能否用迭代来实现 FFT。
我们可以观察一下每次递归，都是哪些元素被分到了一组。接下来拿 \(0\sim 7\) 来举例。
乍一看好像毫无规律，但是我们把最后分的一组中二进制表示出来：
（以下所说的数字全为二进制）
我们可以发现，\(000\) 处在了第 \(000\) 号位，\(100\) 处在了第 \(001\) 号位，\(010\) 处在了第 \(010\) 号位，\(110\) 处在了第 \(011\) 号位......
结论已经显然易见了，把二进制位翻转一下即为 FFT 时的排列顺序，因此我们可以通过这个操作来增加效率，我们把这个操作称为位逆序置换。
不难发现其中有 \(4\) 条规律：
- 奇数位二进制首位为 \(1\)，偶数位二进制首位为 \(0\)。
- 每两个二进制数除了首位都相同。
- 每两个数除首位外是上一级子问题的解
- 前一半数除末位外是上一级子问题的解
因此我们可以写出下面的递推代码：rev[i] = (rev[i >> 1] >> 1) | (i & 1) << (len - 1)
蝴蝶操作
在代码中，有一段操作：
a[k] = a1[k] + w * a2[k];
a[k + lim / 2] = a1[k] - w * a2[k];
其中都包含了共同的一项 w * a2[k]，我们可以将这一项单独存储，减小常数。
同时我们也可以将 a1，a2 数组省去，直接在原数组上进行操作。
代码实现：
#define LOCAL
#include <bits/stdc++.h>
using namespace std;
const double PI = acos(-1);
int lim, len;
int rev[N];
void FFT(Complex a[], int opt)
{
for(int i = 0; i < lim; i ++ )
if(i < rev[i])
swap(a[i], a[rev[i]]);
int up = log2(lim);
for(int dep = 1; dep <= up; dep ++ )
{
int m = 1 << dep;
Complex wn = Complex(cos(2.0 * PI / m), opt * sin(2.0 * PI / m));
for(int k = 0; k < lim; k += m)
{
Complex w = Complex(1, 0);
for(int j = 0; j < m / 2; j ++ )
{
Complex t = w * a[k + j + m / 2];
Complex u = a[k + j];
a[k + j] = u + t;
a[k + j + m / 2] = u - t;
w = w * wn;
}
}
}
if(opt == -1)
{
for(int i = 0; i < lim; i ++ ) a[i].a /= lim;
}
}
int n, m;
int main()
{
n = read(), m = read();
for(int i = 0; i <= n; i ++ ) scanf("%lf", &F[i].a);
for(int i = 0; i <= m; i ++ ) scanf("%lf", &G[i].a);
lim = 1;
while(lim <= n + m) lim <<= 1, len ++;
for(int i = 0; i < lim; i ++ )
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len - 1));
FFT(F, 1), FFT(G, 1);
for(int i = 0; i <= lim; i ++ ) F[i] = F[i] * G[i];
FFT(F, -1);
for(int i = 0; i <= n + m; i ++ ) cout << (int)(F[i].a + 0.5) << ' ';
return 0;
}
浙公网安备 33010602011771号