Appearance
快速傅里叶变换能在
多项式的表示法
系数表示法
对于
系数表示法用
点值表示法
两点确定
点值表示法用至少
多项式乘法
离散傅里叶变换(Discrete Fourier Transform,DFT)是将多项式从 系数表示法 转换为 点值表示法 的算法,IDFT 是 DFT 的逆过程。
运用 DFT 和 IDFT 计算多项式乘法
- 对
和 使用 DFT,得到两个点集: - 计算出
的点值表示法: - 使用 IDFT 将其转化为系数表示法。
其中
一般的 DFT 和 IDFT 时间复杂度高达
让我们从一个例子入手,先熟悉一下
普通的算法
现在计算
这个乘法的结果肯定是一个
第一步:计算
第二步:计算
第三步:转化为系数表示法。这一步的方法有很多,可以使用拉格朗日插值法等。总之最后算出来的结果是
当
上述流程存在一个比较现实的问题:在实际的应用场景中,多项式的次数往往很大。如果我们随意地取
所以
但是只有这三个数是远远不够的。去哪里找其它的数呢?
这样的数,数学家们在复数域中找到了无穷多个。
单位复根
形如
被称为实部; 被称为虚部; 是虚数单位, ; 是这个复数的模。
在复平面上,
表示的是复平面内的横坐标; 表示的是复平面内的纵坐标; - 表示实数
的点都在 轴上,所以 轴又称为「实轴」; - 表示纯虚数
的点都在 轴上,所以 轴又称为「虚轴」。
如图 1,在复平面上画一个半径为
如果把圆周角
中学课本告诉了我们复数乘法的规律:幅角相加模相乘。我们知道
根据上述规律,不难发现,
![]() | ![]() |
---|---|
图 1 | 图 2 |
单位复根还具有如下优异的性质:
- 周期性:
; - 消去性:
; - 对称性:
。
证明并不困难。一方面可以用上文中
也能能轻易得证。
单位复根 恰好可以完美地解决 DFT 中取点的问题。代入
但是,单位复根仅解决了精度问题。目前为止,时间复杂度仍然是
多项式分治
对一个
进行如下变换(假定
- 将偶数项留在前面,将奇数项移到后面:
- 对后一半提取公因式
:
设
则
注意这里我们将
将
将
我们发现,
对每个
再根据之前推出的公式
得出
进而得出我们所需要的
设
以下是
可以画出
容易看出,以上分治算法每次都能递归地将规模为
之所以在一开始假定
例如对于
以上就是 FFT 加速算法的全部内容。不过我们只是用它加速了 DFT 的过程。多项式乘法的最后一步,也就是 IDFT,我们仍然没有提及。实际上 IDFT 也可以用同样的算法进行加速。
IDFT
前文所论述的 DFT 算法实际上是在解以下方程组:
令
DFT 的数学本质就是已知
为了实现 IDFT,我们可以很自然地对上式做出以下变换:
那怎么求中间的这个逆矩阵呢?我们一眼就可以盯真出这是个特殊范德蒙矩阵的逆矩阵。
在范德蒙矩阵中,当
再结合 单位复根 的性质就可以得出
再将这个结果代入进去,并还原成方程组:
这个方程组和原先的方程组极为相似。这意味着我们只需调整 DFT 代码中某些参数的正负号,并且将最终的结果除以
模板
cpp
#include <bits/stdc++.h>
using namespace std;
typedef complex<double> Comp;
const double PI = acos(-1);
vector<Comp> DFT(vector<Comp> a, bool invert) {
int n = a.size();
if (n == 1) return a;
vector<Comp> a0(n / 2), a1(n / 2);
for (int i = 0; 2 * i < n; i ++) {
a0[i] = a[2*i];
a1[i] = a[2*i + 1];
}
vector<Comp> y0 = DFT(a0, invert);
vector<Comp> y1 = DFT(a1, invert);
vector<Comp> y(n);
double angle = 2 * PI / n * (invert ? -1 : 1);
Comp w(1), wn(cos(angle), sin(angle));
for (int i = 0; i < n / 2; i++) {
y[i] = y0[i] + w * y1[i];
y[i + n/2] = y0[i] - w * y1[i];
if (invert) {
y[i] /= 2;
y[i + n/2] /= 2;
}
w *= wn;
}
return y;
}
vector<Comp> multiply(vector<Comp> A, vector<Comp> B) {
int n = 1;
while (n < A.size() + B.size())
n *= 2;
A.resize(n);
B.resize(n);
vector<Comp> yA = DFT(A, false);
vector<Comp> yB = DFT(B, false);
vector<Comp> yC(n);
for (int i = 0; i < n; i ++)
yC[i] = yA[i] * yB[i];
vector<Comp> C = DFT(yC, true);
for (int i = 0; i < n; i ++)
C[i] = round(C[i].real());
while (C.size() && ! C.back().real())
C.pop_back();
return C;
}
int main() {
vector<Comp> A = {1, 2, 3}; // Represents the polynomial 1 + 2x + 3x^2
vector<Comp> B = {4, 5}; // Represents the polynomial 4 + 5x
vector<Comp> C = multiply(A, B);
for (auto i : C)
cout << i.real() << ' ';
cout << endl;
return 0;
}