码迷,mamicode.com
首页 > 其他好文 > 详细

【数学】【多项式】快速傅里叶变换(FFT)

时间:2021-03-01 13:07:46      阅读:0      评论:0      收藏:0      [点我收藏+]

标签:注意   isp   函数   形式   计算   toc   连续   递归   应用   

引入

求两个多项式的卷积

Description

给定两个多项式 \(F\left(x\right), G\left(x\right)\) 的系数表示法,求两个多项式的卷积。

如:

\[F\left(x\right) = 2x + 1 \]

\[G\left(x\right) = x^2 + 2x + 1 \]

得到的结果即为:

\[F\left(x\right)G\left(x\right) = \left(2x + 1\right)\left(x^2 + 2x + 1\right) = 2x^3 + 5x^2 + 4x + 1 \]

一种计算方法是将两个多项式中的每一项两两相乘,最后再将结果相加。

可惜,这样实现的复杂度不够优秀,为 \(\mathcal O\left(n^2\right)\) 级别。

因此我们需要寻找更快的计算方法。

多项式的点值表示和系数表示

我们知道,一个 \(n\) 次多项式,可以被 \(n + 1\) 个横坐标互不相同的点唯一确定。这样用 \(n + 1\) 个点的集合来表示一个多项式的方法被称为点值表示。对应的,用每一项的系数构成的集合表示一个多项式的方法就被称为系数表示

例如方才的多项式 \(G\left(x\right)\),系数表示: \(\{1, 2, 1\}\),点值表示(不唯一):\(\{\left(-1, 0\right), \left(0, 1\right), \left(1, 4\right)\}\)

这样,设 \(n\)\(F\left(x\right)G\left(x\right)\) 的次数,下同,也就意味着我们可以通过找 \(n + 1\)\(x\) 坐标,求出 \(F\left(x\right)\)\(G\left(x\right)\) 在这,\(n + 1\) 个位置上的函数值,再将这 \(n - 1\) 个位置上的函数值分别相乘,就得到了 \(F\left(x\right)G\left(x\right)\) 的一个点值表示。

计算的复杂度是 \(\mathcal O\left(n\right)\) 级别,非常优秀。

求值和插值

实际应用上,我们很少通过点值表示来表示一个多项式。

因此,若想快速计算两个多项式的卷积,需要将两个多项式由系数表示转换为点值表示,然后在点值表示下计算两个多项式的卷积,再将点值表示转换为系数表示。

计算卷积的过程已经解决,现在需要的是实现将系数表示转换为点值表示,和将点值表示转换为系数表示的方法。这两个过程分别被称为求值插值

考虑暴力计算。找到 \(n + 1\) 个两两不同的 \(x\) 值,然后分别代入两个多项式里去,求出对应的点值,进行完乘法后,再用待定系数法,根据得到的点列出 \(n + 1\) 个方程,然后消元,或者使用插值法。复杂度?

消元复杂度:\(\mathcal O\left(n^3\right)\);插值复杂度:单个 \(\mathcal O\left(n^2\right)\),在 \(x\) 取值连续时可以优化到 \(\mathcal O\left(n\right)\),由于 \(x\) 值可以任取,可以认为其复杂度为 \(\mathcal O\left(n^2\right)\)。没有起到优化作用。

快速傅里叶变换(FFT)

一种思路

考虑一个多项式 \(F\left(x\right) = x^2\)。发现因为它的函数图像具有对称性,因此只需要取 \(\dfrac{n}{2}\)\(x\) 值,就可以得到 \(n\) 个点值。

那对于 \(F\left(x\right) = x^3\) 呢?显然还是有的。

那对于一个任意多项式 \(P\left(x\right)\) 呢?

(接下来我们默认多项式的项数为 \(2^n\) 形式,即便不是,也可以简单地将不足的项的系数视为 \(0\)。)

考虑将多项式 \(P\left(x\right)\)

\[P\left(x\right) = \sum_{i = 0} ^ n a_i x^i \]

将其按奇偶次项分成两个多项式 \(P_e\left(x\right), P_o\left(x\right)\)

\[P_e\left(x\right) = \sum_{i = 0} ^ {\frac{n}{2}} a_{2i} \left(x^2\right)^i, P_o\left(x\right) = \sum_{i = 0} ^ {\frac{n}{2}} a_{2i + 1} x^{2i + 1} = x\sum_{i = 0} ^ {\frac{n}{2}} a_{2i + 1} \left(x^2\right)^i \]

则显然有:

\[P\left(x\right) = \sum_{i = 0} ^ {\frac{n}{2}} a_{2i} \left(x^2\right)^i + x\sum_{i = 0} ^ {\frac{n}{2}} a_{2i + 1} \left(x^2\right)^i = P_e\left(x\right) + P_o\left(x\right) \]

\[P\left(-x\right) = \sum_{i = 0} ^ {\frac{n}{2}} a_{2i} \left(x^2\right)^i - x\sum_{i = 0} ^ {\frac{n}{2}} a_{2i + 1} \left(x^2\right)^i = P_e\left(x\right) - P_o\left(x\right) \]

看,发现了什么?我们只需要求出 \(P_e\left(x\right)\)\(P_o\left(x\right)\),就可以得到 \(P\left(x\right)\)\(P\left(-x\right)\) 两处的多项式的点值!

\(P_e\left(x\right), P_o\left(x\right)\) 又是两个关于 \(x^2\)\(\frac{n}{2}\) 次多项式,这似乎意味着我们可以在 \(\mathcal O\left(n \log n\right)\) 的时间复杂度级别下,递归地来实现这个算法!

但是,注意到这里面有什么问题吗?因为这个算法能够实现,是基于求 \(P\left(x\right)\) 的点值时 \(x\) 的取值能够正负配对的假设。而 \(x^2, x^4\) 等偶次项,不是正负配对的。这样递归就无法进行。

解决方案

都进行到这里了,有没有解决方案呢?答案是肯定的。

在实数域内,显然是不能使得 \(x\) 的取值总是能够两两配对的了。因此考虑把 \(x\) 的取值扩展到复数域。

我们知道, \(\mathbb i^2 = -1\)。因此若 \(x\) 的取值集合是 \(\{1, -1, \mathbb i, -\mathbb i\}\) 的话,\(x^2\) 就也满足两两配对了。

而这四个数又分别是什么呢?没错,就是 \(4\) 次单位根。

这便是 FFT 的基本思想。

快速傅里叶变换(DFT)

通过系数表示求点值表示的方式称为 DFT。

\(n\) 次单位根(\(\omega_n\))可以通过如下的公式求得,其中 \(\alpha\) 是幅角,采用弧度制,证明自行搜索复数相关内容:

\[\omega_n = \cos\left(\alpha\right) + i\sin\left(\alpha\right) \]

或者如下的公式求得,其中 \(e\) 是自然对数的底数:

\[\omega_n = e^{\frac{2\pi i}{n}} \]

根据单位根的性质有:

\[\omega_{2n}^{2k} = \omega_n^k, \omega_{n}^{k + \frac{n}{2}} = - \omega_n ^ k \]

这样,想要求 \(P\left(x\right)\)\(\{1, \omega_n, \omega_n^2, \cdots \omega_n^{\left(n - 1\right)}\}\) 时的点值,转而求 \(P_e\left(x\right)\)\(P_o\left(x\right)\) 分别在 \(\{1, \omega_{\frac{n}{2}}, \omega_{\frac{n}{2}} ^ 2, \cdots, \omega_{\frac{n}{2}} ^ {\frac{n}{2} - 1}\}\) 时的点值,这样递归下去,在要求的多项式只有一项的时候返回系数,合并的时候把得到的 \(P_o\left(x\right)\) 的结果乘上一个当前的单位根,然后再分别一加一减计算出 \(P\left(x\right)\)\(P\left(-x\right)\)

伪代码:

DFT(P{p_0, p_1, p_2, ..., p_n-1}):
  if n = 1
    返回 P 
  定义两个求值点子序列 Pe{p_0, p_2,...} Po{p_1, p_3,...}
  DFT(Pe) DFT(Po)
  计算 n 次单位根 wn
  枚举 i
    将 Pe_i + wn ^ i * Po_i 和 Pe_i - wn ^ i * Po_i 填入结果序列对应的位置

快速傅里叶逆变换(IDFT)

求值的问题解决了,我们现在来解决插值。

不难想到,插值应该和求值是基本相同的问题。

考虑 DFT 前的系数矩阵,让它乘上 DFT 矩阵后,变为点值矩阵,则 DFT 矩阵即运算为:

\[\begin{bmatrix}v_0 \\ v_1 \\ v_2 \\ v_3 \\ \vdots \\ v_{n-1} \end{bmatrix} = \begin{bmatrix}1 & 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_n^1 & \omega_n^2 & \omega_n^3 & \cdots & \omega_n^{n-1} \\ 1 & \omega_n^2 & \omega_n^4 & \omega_n^6 & \cdots & \omega_n^{2(n-1)} \\ 1 & \omega_n^3 & \omega_n^6 & \omega_n^9 & \cdots & \omega_n^{3(n-1)} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \omega_n^{3(n-1)} & \cdots & \omega_n^{(n-1)^2} \end{bmatrix} \times \begin{bmatrix} c_0 \\ c_1 \\ c_2 \\ c_3 \\ \vdots \\ c_{n-1} \end{bmatrix} \]

奇妙的,若找到一个矩阵,使得它乘上点值矩阵能够得到系数矩阵,则这个矩阵就是 IDFT 矩阵。

不难想到 IDFT 矩阵是 DFT 矩阵的逆矩阵。

那么 IDFT 矩阵是什么样子的呢?非常奇妙:

\[\begin{bmatrix} 1 & 1 & 1 & 1 & \cdots & 1 \1 & \omega_n^{-1} & \omega_n^{-2} & \omega_n^{-3} & \cdots & \omega_n^{-(n-1)} \1 & \omega_n^{-2} & \omega_n^{-4} & \omega_n^{-6} & \cdots & \omega_n^{-2(n-1)} \1 & \omega_n^{-3} & \omega_n^{-6} & \omega_n^{-9} & \cdots & \omega_n^{-3(n-1)} \\vdots & \vdots & \vdots & \vdots & \ddots & \vdots \1 & \omega_n^{-(n-1)} & \omega_n^{-2(n-1)} & \omega_n^{-3(n-1)} & \cdots & \omega_n^{-(n-1)^2} \end{bmatrix}\]

也就是说,要实现 IDFT,只需要将原来 DFT 的代码中,改变函数名称,然后将单位根改成单位根的 \(-1\) 次方即可!

【数学】【多项式】快速傅里叶变换(FFT)

标签:注意   isp   函数   形式   计算   toc   连续   递归   应用   

原文地址:https://www.cnblogs.com/zimujun/p/14456655.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!