取自 浅谈 FFT (1){:target=”_blank”}
FFT 快速傅里叶变换
理论知识
回归最开始讨论的内容,多项式
同上,默认 $n$ 为 $2$ 的正整数次幂,不足则在多项式后面补 $0$
设一个 $n$ 次多项式 $A(x) = a _ 0 + a _ 1 x + a _ 2 x ^ 2 \cdots + a _ {n - 2} x ^ {n - 2} + a _ {n - 1} x ^ {n - 1}$
将其每项按所在位置的奇偶性分成两部分,得
$A(x) = (a _ 0 + a _ 2 x ^ 2 + \cdots + a _ {n - 2} x ^ {n - 2}) + (a _ 1 x + a _ 3 x ^ 3 + \cdots + a _ {n - 1} x ^ {n - 1}) $
设
$A _ 1(x) = a _ 0 + a _ 2 x ^ 2 + \cdots + a _ {n - 2} x ^ {\frac n 2 - 1}$
$A _ 2(x) = a _ 1 + a _ 3 x ^ 1 + \cdots + a _ {n - 1} x ^ {\frac n 2 - 1}$
那么有
$ \begin{aligned} A (x) = A _ 1(x ^ 2) + x A _ 2(x ^ 2) \\\ \end{aligned} $
假设 $0 \leq k < \frac n 2$,将 $x = \omega _ n ^ k$ 代入,得
$A(\omega _ n ^ k) = A _ 1(\omega _ n ^ {2 k}) + \omega _ n ^ k A _ 2(\omega _ n ^ {2 k})$
将 $x = \omega _ n ^ {k + \frac n 2}$ 代入,得
$ \begin{aligned} A(\omega _ n ^ {k + \frac n 2}) = &\ A _ 1(\omega _ n ^ {2 k + n}) + \omega _ n ^ {k + \frac n 2} A _ 2(\omega _ n ^ {2 k + n}) \\\ \\\ = &\ A _ 1(\omega _ n ^ {2 k}) + \omega _ n ^ {k + \frac n 2} A _ 2(\omega _ n ^ {2 k}) \quad (根据性质 5) \\\ \\\ = &\ A _ 1(\omega _ n ^ {2 k}) - \omega _ n ^ k A _ 2(\omega _ n ^ {2 k}) \quad (根据性质 4) \\\ \\\ \end{aligned} $
发现了吗?
没发现的话,对比一下再看看?
$ \begin{aligned} A(\omega _ n ^ k) = &\ A _ 1(\omega _ n ^ {2 k}) + \omega _ n ^ k A _ 2(\omega _ n ^ {2 k}) \\\ A(\omega _ n ^ {k + \frac n 2}) = &\ A _ 1(\omega _ n ^ {2 k}) - \omega _ n ^ {k} A _ 2(\omega _ n ^ {2 k}) \\\ \end{aligned} $
没错!这两个式子只有一个符号不同!
那也就是说,我们在求 $A(\omega _ n ^ k)$ 的时候,可以 $O(1)$ 得到第二个式子的值!
而第一个式子,在 $k$ 取遍 $[0 \sim \frac n 2 - 1]$ 时,第二个式子正好取遍 $[\frac k 2 \sim k]$!
那么我们就只需要将 $[0 \sim \frac n 2 - 1]$ 的所有值代入 $A$ 求值即可,
原问题缩小了一半!!
不仅如此,子问题还是满足原问题的性质的,所以我们就可以递归分治去做!!
递归边界:整个式子只剩一项,直接返回即可
时间复杂度分析
在每层递归中,将原函数的奇数项和偶数项分别处理,
所以最多只会递归 $\log n$ 层
每层递归中
- 要将原函数的奇数项和偶数项分离开,时间复杂度为 $O(n)$
- 用得到的值,求出 $A(\omega _ n ^ k)$ 和 $A(\omega _ n ^ {k + \frac n 2})$ 的值 $(0 \leq k < \frac n 2)$,时间复杂度 $O(n)$
所以总的时间复杂度为 $\log n \times O(n) = O(n \log n)$
C ++ 代码实现
代码中要用到虚数,可以直接使用 STL
的 complex
但不开 O2
可能会被卡
这里先给出手写 complex
的板子
由于在 $FFT$ 中暂时不需要除法,所以就先不写除法了
typedef struct complex
{
// a 为实部,b 为虚部
double a, b;
complex(double x = 0, double y = 0)
{
a = x, b = y;
}
// 重载加法
complex operator + (const complex &t) const
{
return complex(a + t.a, b + t.b);
}
// 重载减法
complex operator - (const complex &t) const
{
return complex(a - t.a, b - t.b);
}
// 重载乘法,上面 2.2.5 的乘法法则中有推导过程
complex operator * (const complex &t) const
{
return complex(a * t.a - b * t.b, a * t.b + b * t.a);
}
} cp;
下面是用递归FFT将系数表示法转换成点值表示法的代码
这段代码跑完后,a[i]
存的是 $A(\omega _ n ^ i)$ 的值
// 将 n 次多项式 a 从系数表示转为点值表示
void FFT(cp a[], int n)
{
// 如果 a 中只剩一项,直接返回
if (n == 1) return;
cp t[N]; // 临时数组
int m = n >> 1; // 分成相等长度的两个多项式处理
for (int i = 0; i < m; i ++ )
t[i] = a[i << 1], // t[0 ~ m) 存偶数项
t[i + m] = a[i << 1 | 1]; // t[m ~ n) 存奇数项
FFT(t, m); // 分治处理偶数项
FFT(t + m, m); // 分治处理奇数项
// 这两个递归做完后,
// t[0 ~ m) 中,t[i] 存的是 A_1(omega_n ^ 2i) 的值 (0 <= k < m = n / 2)
// t[m ~ n) 中,t[i] 存的是 A_2(omega_n ^ 2i) 的值 (0 <= k < m = n / 2)
// 合并到 a 中
for (int i = 0, j = m; i < m; i ++ , j ++ )
{
// 计算 omega_n ^ i 的值,存入 x 中
cp x(cos(2 * pi * i / n), sin(2 * pi * i / n));
// 根据上面两个式子,分别求出 k 取 1 ~ n/2 时,
// A(omega_n ^ k) 和 A(omega_n ^ (k + n/2)) 的值
a[i] = t[i] + x * t[j];
a[j] = t[i] - x * t[j];
}
}
IFFT 快速傅里叶逆变换
别忘了,我们现在只解决了将系数表示法转换成点值表示法的问题
我们将两个多项式转换到点痣表示法后,再做乘积,得到的是点值表示法下的乘积
但我们平时用系数表示的多项式,对我们而言,点值表示的多项式没有意义
现在我们要考虑如何快速地将点值表示法转换成系数表示法,这个过程叫做傅里叶逆变换
这里先给出结论,再给出推导过程。因为这个结论的确很难想出来
理论知识
结论:
设函数
$ \begin{aligned} A(x) = a _ 0 + a _ 1 x + a _ 2 x ^ 2 + \cdots + a _ {n - 2} x ^ {n - 2} + a _ {n - 1} x _ {n - 1} \end{aligned} $
那么有
$ a _ i = \displaystyle \frac { \begin{aligned} \sum _ {j = 0} ^ {n - 1} A(\omega _ n ^ j) \omega _ n ^ {- i j} \end{aligned} } n $
证明:
这个证明比较高能,请做好心理准备
推式子,先将分子那项提出来
$ \begin{aligned} & \sum _ {j = 0} ^ {n - 1} A(\omega _ n ^ j) \omega _ n ^ {- i j} \\\ \\\ & = \sum _ {j = 0} ^ {n - 1} \omega _ n ^ {- i j} \sum _ {k = 0} ^ {n - 1} a _ k (\omega _ n ^ j) ^ k \qquad (将 A(\omega _ n ^ j) 拆开) \\\ \\\ & = \sum _ {j = 0} ^ {n - 1} \sum _ {k = 0} ^ {n - 1} a _ k \omega _ n ^ {j k} \omega _ n ^ {- i j} \qquad (将 \omega _ n ^ {- i j} 乘到里面去,(\omega _ n ^ j) ^ k 指数展开) \\\ \\\ & = \sum _ {j = 0} ^ {n - 1} \sum _ {k = 0} ^ {n - 1} a _ k (\omega _ n ^ {k - i}) ^ j \qquad (将 \omega _ n ^ {j k} \omega _ n ^ {- i j} 合并) \\\ \\\ & = \begin{matrix} a _ 0 (\omega _ n ^ {- i}) ^ 0 + a _ 1 (\omega _ n ^ {1 - i}) ^ 0 + a _ 2 (\omega _ n ^ {2 - i}) ^ 0 + \cdots + a _ {n - 2} (\omega _ n ^ {n - 2 - i}) ^ 0 + a _ {n - 1} (\omega _ n ^ {n - 1 - i}) ^ 0 + \\\ a _ 0 (\omega _ n ^ {- i}) ^ 1 + a _ 1 (\omega _ n ^ {1 - i}) ^ 1 + a _ 2 (\omega _ n ^ {2 - i}) ^ 1 + \cdots + a _ {n - 2} (\omega _ n ^ {n - 2 - i}) ^ 1 + a _ {n - 1} (\omega _ n ^ {n - 1 - i}) ^ 1 + \\\ a _ 0 (\omega _ n ^ {- i}) ^ 2 + a _ 1 (\omega _ n ^ {1 - i}) ^ 2 + a _ 2 (\omega _ n ^ {2 - i}) ^ 2 + \cdots + a _ {n - 2} (\omega _ n ^ {n - 2 - i}) ^ 2 + a _ {n - 1} (\omega _ n ^ {n - 1 - i}) ^ 2 + \\\ + \\\ \vdots \\\ + \\\ a _ 0 (\omega _ n ^ {- i}) ^ {n - 2} + a _ 1 (\omega _ n ^ {1 - i}) ^ {n - 2} + a _ 2 (\omega _ n ^ {2 - i}) ^ {n - 2} + \cdots + a _ {n - 1} (\omega _ n ^ {n - 2 - i}) ^ {n - 2} + a _ {n - 1} (\omega _ n ^ {n - 1 - i}) ^ {n - 2} \\\ a _ 0 (\omega _ n ^ {- i}) ^ {n - 1} + a _ 1 (\omega _ n ^ {1 - i}) ^ {n - 1} + a _ 2 (\omega _ n ^ {2 - i}) ^ {n - 1} + \cdots + a _ {n - 2} (\omega _ n ^ {n - 2 - i}) ^ {n - 1} + a _ {n - 1} (\omega _ n ^ {n - 1 - i}) ^ {n - 1} \\\ \end{matrix} (展开) \\\ \\\ & = \begin{matrix} a _ 0 (\omega _ n ^ {- i}) ^ 0 + a _ 0 (\omega _ n ^ {- i}) ^ 1 + a _ 0 (\omega _ n ^ {- i}) ^ 2 + \cdots + a _ 0 (\omega _ n ^ {- i}) ^ {n - 2} + a _ 0 (\omega _ n ^ {- i}) ^ {n - 1} \\\ a _ 1 (\omega _ n ^ {1 - i}) ^ 0 + a _ 1 (\omega _ n ^ {1 - i}) ^ 1 + a _ 1 (\omega _ n ^ {1 - i}) ^ 2 + \cdots + a _ 1 (\omega _ n ^ {1 - i}) ^ {n - 2} + a _ 1 (\omega _ n ^ {1 - i}) ^ {n - 1} \\\ a _ 2 (\omega _ n ^ {2 - i}) ^ 0 + a _ 2 (\omega _ n ^ {2 - i}) ^ 1 + a _ 2 (\omega _ n ^ {2 - i}) ^ 2 + \cdots + a _ 2 (\omega _ n ^ {2 - i}) ^ {n - 2} + a _ 2 (\omega _ n ^ {2 - i}) ^ {n - 1} \\\ + \\\ \vdots \\\ + \\\ a _ {n - 2} (\omega _ n ^ {n - 2 - i}) ^ 0 + a _ {n - 2} (\omega _ n ^ {n - 2 - i}) ^ 1 + a _ {n - 2} (\omega _ n ^ {n - 2 - i}) ^ 2 + \cdots + a _ {n - 2} (\omega _ n ^ {n - 2 - i}) ^ {n - i - 1} + a _ {n - 2} (\omega _ n ^ {n - 2 - i}) ^ {n - i} + \\\ a _ {n - 1} (\omega _ n ^ {n - 1 - i}) ^ 0 + a _ {n - 1} (\omega _ n ^ {n - 1 - i}) ^ 1 + a _ {n - 1} (\omega _ n ^ {n - 1 - i}) ^ 2 + \cdots + a _ {n - 1} (\omega _ n ^ {n - 1 - i}) ^ {n - i - 1} + a _ {n - 1} (\omega _ n ^ {n - 1 - i}) ^ {n - i} \\\ \end{matrix} (转置) \\\ \\\ & = \begin{matrix} a _ 0 \sum _ {j = 0} ^ {n - 1} (\omega _ n ^ {- i}) ^ j + \\\ a _ 1 \sum _ {j = 0} ^ {n - 1} (\omega _ n ^ {1 - i}) ^ j + \\\ a _ 2 \sum _ {j = 0} ^ {n - 1} (\omega _ n ^ {2 - i}) ^ j + \\\ + \\\ \vdots \\\ + \\\ a _ {n - 2} \sum _ {j = 0} ^ {n - 1} (\omega _ n ^ {n - 2 - i}) ^ j \\\ a _ {n - 1} \sum _ {j = 0} ^ {n = 1} (\omega _ n ^ {n - 1 - i}) ^ j \\\ \end{matrix} (将每行提取公因式,并写成求和的形式) \\\ \\\ & = \sum _ {k = 0} ^ {n - 1} a _ k (\sum _ {j = 0} ^ {n - 1} (\omega _ n ^ {k - i}) ^ j) \qquad (将整个列写成求和的形式) \\\ \end{aligned} $
观察 $ \begin{aligned} \sum _ {j = 0} ^ {n - 1} (\omega _ n ^ {k - i}) ^ j \\\ \end{aligned} $ 这项,
设 $\begin{aligned} S = \sum _ {j = 0} ^ {n - 1} (\omega _ n ^ {k - i}) ^ j = 1 + \omega _ n ^ {k - i} + (\omega _ n ^ {k - i}) ^ 2 + \cdots + (\omega _ n ^ {k - i}) ^ {n - 2} + (\omega _ n ^ {k - i}) ^ {n - 1} \\\ \end{aligned} $
那么有 $\begin{aligned} \omega _ n ^ {k - i} S = &\ \omega _ n ^ {k - i} + (\omega _ n ^ {k - i}) ^ 2 + (\omega _ n ^ {k - i}) ^ 3 \cdots + \cdots (\omega _ n ^ {k - i}) ^ {n - 1} + (\omega _ n ^ {k - i}) ^ n \\\ = &\ S - 1 + (\omega _ n ^ {k - i}) ^ n \\\ \end{aligned}$
得 $\omega _ n ^ {k - i} S = S + (\omega _ n ^ {k - i}) ^ n - 1$
$(\omega _ n ^ {k - i} - 1) S = (\omega _ n ^ {k - i}) ^ n - 1$
当 $\omega _ n ^ {k - i} = 1$ 时,即 $k \equiv i \ (\text{mod}\ n)$ 时,
$S = 1 + 1 + 1 + \cdots + 1 = n$
当 $\omega _ n ^ {k - i} \neq 1$ 时,即 $k \not \equiv i \ (\text{mod}\ n)$ 时,可以将两边同时 $\div (\omega _ n ^ {k - i} - 1)$,得
$ \begin{aligned} S = &\ \frac {(\omega _ n ^ {k - i}) ^ n - 1} {\omega _ n ^ {k - i} - 1} \\\ \\\ = &\ \frac {(\omega _ n ^ n) ^ {k - i} - 1} {\omega _ n ^ {k - i} - 1} \\\ \\\ = &\ \frac {1 ^ {k - i} - 1} {\omega _ n ^ {k - i} - 1} \\\ \end{aligned} $
注意到此时分子为零,且分母不为零,故 $S = 0$
所以有 $\begin{aligned} \sum _ {j = 0} ^ {n - 1} (\omega _ n ^ {k - i}) ^ j = \left\\{ \begin{aligned} & n & \text{if}\ k \equiv i\ (\text{mod}\ n) \\\ & 0 & \text{if}\ k \not \equiv i\ (\text{mod}\ n) \\\ \end{aligned} \right. \end{aligned} $
那么在 $\begin{aligned} \sum _ {k = 0} ^ {n - 1} a _ k (\sum _ {j = 0} ^ {n - 1} (\omega _ n ^ {k - i}) ^ j) \end{aligned}$ 中,对于后面求和那项的值
当 $k = i$ 时,值为 $n$
当 $k \neq i$ 时,值为 $0$
所以原式 $
\begin{aligned}
\\\
\displaystyle \frac
{
\begin{aligned}
\sum _ {j = 0} ^ {n - 1} A(\omega _ n ^ j) \omega _ n ^ {- i j}
\end{aligned}
} n = &\
\frac {
\begin{aligned}
\sum _ {k = 0} ^ {n - 1} a _ k (\sum _ {j = 0} ^ {n - 1} (\omega _ n ^ {k - i}) ^ j)
\end{aligned}
} n \\\
= &\ \frac {n \times a _ i} n \\\
= &\ a _ i \\\
\end{aligned}
$
证毕
时间复杂度分析
由于我们是已经知道了所有 $A(\omega _ n ^ j)$ 的值的,
所以我们在求解 $\begin{aligned} \sum _ {j = 0} ^ {n - 1} A(\omega _ n ^ j) \omega _ n ^ {- i j} \end{aligned}$ 的时候呢,
那么我们就可以将这个式子中的 $A(\omega _ n ^ j)$ 看做是整个多项式的系数,
而将 $\omega _ n ^ {-i}$ 看做是整个多项式的自变量,
那么我们就可以通过快速傅里叶变换,做这个快速傅里叶逆变换
时间复杂度为 $O(n \log n)$
C ++ 代码实现
void FFT(cp a[], int n)
{
if (n == 1) return;
cp t[N];
int m = n >> 1;
for (int i = 0; i < m; i ++ )
t[i] = a[i << 1],
t[i + m] = a[i << 1 | 1];
FFT(t, m);
FFT(t + m, m);
for (int i = 0, j = m; i < m; i ++ , j ++ )
{
// 这里只需要将后面的 sin(2 pi i / n) 改为 -sin(2 * pi * i / n) 即可
// 这么改的原因,是因为我们在该函数中,代入的变量变成了 omega_n ^ -i
// 而 omega_n ^ -i = cos(2 * pi * -i / n) + sin(2 * pi * -i / n)
// = cos(2 * pi * i / n) - sin(2 * pi * i / n) (诱导公式三)
cp x(cos(2 * pi * i / n), -sin(2 * pi * i / n));
a[i] = t[i] + x * t[j];
a[j] = t[i] - x * t[j];
}
}
多项式乘法的递归实现
C ++ 代码
#include <cmath>
#include <cstdio>
#include <cstring>
typedef struct complex
{
double a, b;
complex(double x = 0, double y = 0)
{
a = x, b = y;
}
complex operator + (const complex &t) const
{
return complex(a + t.a, b + t.b);
}
complex operator - (const complex &t) const
{
return complex(a - t.a, b - t.b);
}
complex operator * (const complex &t) const
{
return complex(a * t.a - b * t.b, a * t.b + b * t.a);
}
} cp;
const int N = 400005;
const double pi = 2 * acos(-1);
int n, m;
cp a[N], b[N];
void FFT(cp a[], int n, int f)
{
if (n == 1) return;
static cp t[N];
int m = n >> 1;
for (int i = 0; i < m; i ++ )
t[i] = a[i << 1],
t[i + m] = a[i << 1 | 1];
// 将 t 复制到 a 中
// 这里偷个懒,就直接用 memcpy 了
// 一个 double 占 8 个字节
// complex 有两个 double,占了 16 个字节
// 所以这里要将 n 乘 16
memcpy(a, t, n << 4);
// 然后对 a 做剩下的 FFT
FFT(a, m, f);
FFT(a + m, m, f);
for (int i = 0, j = m; i < m; i ++ , j ++ )
{
cp x(cos(pi * i / n), f * sin(pi * i / n));
t[i] = a[i] + x * a[j]; // 这里要改成从 a 到 t 做
t[j] = a[i] - x * a[j];
}
// 然后再将 t 复制回 a
memcpy(a, t, n << 4);
}
int main()
{
scanf("%d %d", &n, &m);
for (int i = 0; i <= n; i ++ )
{
double x;
scanf("%lf", &x);
a[i] = {x, 0};
}
for (int i = 0; i <= m; i ++ )
{
double x;
scanf("%lf", &x);
b[i] = {x, 0};
}
int k = 1;
while(k <= n + m) k *= 2;
FFT(a, k, 1);
FFT(b, k, 1);
for (int i = 0; i < k; i ++ ) a[i] = a[i] * b[i];
FFT(a, k, -1);
for (int i = 0; i <= n + m; i ++ )
printf("%d ", (int)(a[i].a / k + 0.5));
return 0;
}
妙!
抽风大佬tql
%%%%,谢谢抽风大佬整理板书
请问n不为2的正整数次幂怎么办呢?一直很疑惑QAQ
突然明白,好像在后面补系数为0的项就行了(
证明让我直接劝退…
tql
%%%