实序列离散傅里叶变换的快速算法

January 21, 2018January 25, 2018 miskcoo 2 Comments

$\newcommand{\dft}[1]{\text{DFT}\left[#1\right]}$ 在从多项式乘法到快速傅里叶变换中我们介绍了快速傅立叶变换的算法。由于在计算多项式乘法或者大整数乘法时，我们所使用的都是实数，而快速傅立叶变换并不仅适用于实序列，对于复序列也是可以计算的，这让我们考虑是否可以利用虚部来保存额外的信息而减少计算的次数。

事实上，对于实序列 $x(n)$ ，其DFT是共轭对称的，也就是 $X^*(n) = X(N - n)$ ，如果将两个实数序列 $x(n)$ 和 $y(n)$ 合并称为新的 $z(n) = x(n) + i\cdot y(n)$ ，对其进行FFT，我们可以从中复原出分别对两个序列进行FFT的结果。也就是，可以通过两次FFT就进行多项式乘法的计算（原先是需要三次的）。

沿用其中的记号，先回忆一下离散傅里叶变换的公式。

考虑长度为 $N$ 的序列 $x(n)$ 及其 $N$ 点离散傅里叶变换 $X(n) := \dft{x(n)}$

\begin{equation} X(n) = \sum_{k=0}^{N - 1} x(k)\cdot\omega_N^{kn},~ n = 0, 1, \cdots, N - 1 \end{equation}

其中 $\omega_N = e^{-2\pi i/N}$ ，即 $N$ 次单位根。

现在来考虑 $X(n)$ 和 $X(N - n)$ 之间的关系，我们用 $X^*(n)$ 或 $\overline{X(n)}$ 来表示 $X(n)$ 的共轭

\begin{equation} X^*(n) = \sum_{k=0}^{N - 1} x^*(k)\cdot\omega_N^{-kn} = \sum_{k=0}^{N - 1} x^*(k)\cdot\omega_N^{k(N-n)} \end{equation}

这是因为 $\overline{\omega_N^{kn}} = e^{2\pi ikn/N} = e^{-2\pi ik(N - n)/N} = \omega_N^{k(N - n)}$

那么，对于实序列， $x^*(n) = x(n)$ ，我们可以得到 $X^*(n) = X(N - n)$ 。也就是说，实序列的DFT是共轭对称的。

现在来考虑两个实序列 $x(n)$ 和 $y(n)$ ，其对应的DFT为 $X(n)$ 和 $Y(n)$ ，定义新的序列 $z(n) = x(n) + i\cdot y(n)$ ，其对应DFT为 $Z(n)$ ，那么

\begin{equation} Z(n) = \dft{x(n) + i\cdot y(n)} = X(n) + i\cdot Y(n) \end{equation}

注意这里 $X(n)$ 和 $Y(n)$ 本身并不是一个实数，因此它们也不是 $Z(n)$ 的实部和虚部，我们设 $Z(n)$ 的实部和虚部分别为 $A(n)$ 和 $B(n)$ 。

\begin{equation} \begin{aligned} \overline{X(N - n) + i\cdot Y(N - n)} &= \overline{A(N - n) + i\cdot B(N - n)} \\ \Rightarrow X(n) - i\cdot Y(n) &= A(N - n) - i\cdot B(N - n) \end{aligned} \end{equation}

这样就可以得到

\begin{equation} \begin{aligned} 2X(n) &= [A(n) + A(N - n)] + i[B(n) - B(N - n)] \\ 2Y(n) &= [B(n) + B(N - n)] - i[A(n) - A(N - n)] \end{aligned} \end{equation}

这样，对于两个实序列的DFT可以只利用一次FFT以及一点额外的运算进行计算而不必分别计算两次FFT。

下面来考虑一个长度为 $2N$ 的实序列 $x(n)$ 。我们将其按照下标的奇偶分解为两个长度为 $N$ 的序列

\begin{equation} \begin{aligned} x_1(n) &= x(2n)\\ x_2(n) &= x(2n + 1) \end{aligned} \end{equation}

现在可以利用刚刚的方法通过一次N点的FFT来计算它们分别的DFT序列 $X_1(n)$ 和 $X_2(n)$ ，通过类似于FFT的办法将其合并

\begin{equation} \begin{aligned} X_1(n) + \omega_{2N}^n X_2(n) &= \sum_{k=0}^{N-1}x_1(k)\omega_N^{kn}+\omega_{2N}^n\sum_{k=0}^{N-1}x_2(k)\omega_N^{kn} \\ &= \sum_{k=0}^{N-1}x(2k)\omega_{2N}^{2kn}+\sum_{k=0}^{N-1}x(2k + 1)\omega_{2N}^{(2k+1)n} \\ &= X(n) \end{aligned} \end{equation}

这里只有前 $N$ 项，对于其余项可以利用共轭对称这个性质直接计算。对于大整数乘法，我们也可以利用两个实序列DFT的快速算法进行计算，这样只需要两次FFT即可，下面我提供一个多项式乘法的代码（UOJ#34），FFT的模板是很久以前写的就别看了，可以看看下面处理的部分。

(more…)

特殊多项式在整点上的线性插值方法

August 4, 2015August 19, 2017 miskcoo 2 Comments

首先我来解释一下题目是什么东西，假如给你一个 $k$ 次多项式 $P(x)$ ，并且已知了 $P(0)$ , $P(1)$ , $\cdots$ , $P(k)$ ，要求计算出 $P(n), n \in \mathbb N$ 的值，由于已经知道了 $k + 1$ 个点值，那么要插值出系数是可以在 $\mathcal O(k^2)$ 完成，如果使用 FFT 的话是可以在 $\mathcal O(k \log^3 k)$ 的时间内完成的^[1]，但是对于这个特别的问题，我们有一个 $\mathcal O(k)$ 的算法！这好像是杜教先提出来的，然后还有这篇文章

(more…)

牛顿迭代法在多项式运算的应用

June 8, 2015August 19, 2017 miskcoo 3 Comments

总算是快把 FFT 和生成函数的各种东西补了好多，膜拜策爷的论文和 Picks 的博客 QAQ

这篇文章大概就是说如何用牛顿迭代法来解满足 $G(F(z)) \equiv 0 \pmod {z^n}$ 的 $F(z)$

然后这东西可以比较方便地计算 $\sqrt{A(z)}$ 、 $e^{A(z)}$ ，也就是多项式开根、求指数对数之类鬼畜的东西，在生成函数计数中十分有用

顺便一提，这里说的“多项式”实际上你可以直接理解为生成函数或者形式幂级数

(more…)

多项式的多点求值与快速插值

May 25, 2015August 19, 2017 miskcoo 7 Comments

多项式的多点求值（multi-point evaluation）是给出一个多项式 $A(x)$ ，和 $n$ 个点 $x_0, x_1, \cdots, x_{n-1}$ ，要求求出 $A(x_0), A(x_1), \cdots, A(x_{n-1})$

相反的，多项式的插值（interpolation）是给出 $n+1$ 个点 $(x_0, y_0), (x_1, y_1), \cdots, (x_n, y_n)$ ，求出一个 $n$ 次多项式，使得这些点都在这个多项式上

这两个问题实际上是在多项式的点值表示（point-value representation）和系数表示（coefficient representation）之间转换的方法，在快速傅里叶变换中由于带入值的特殊性质，可以在 $\mathcal O(n\log n)$ 的之间内将两种东西互相转换，但是，如果是任意给定点要求求值或者插值，就没有比较好的性质可以利用，但是仍然有比较快速的方法来计算它们

(more…)

多项式求逆元

May 10, 2015August 19, 2017 miskcoo 11 Comments

概述

多项式求逆元是多项式除法和多项式取模的必要过程，在竞赛中，有时候多项式求逆元的应用比多项式的除法还要多，用快速傅里叶变换以及倍增算法可以做到在 $\mathcal O(n\log n)$ 的时间内求出一个多项式的逆元

基本概念

在介绍多项式的逆元之前，先说明一些概念：多项式的度、多项式的逆元、多项式的除法和取余

对于一个多项式 $A(x)$ ，称其最高项的次数为这个多项式的度（degree），记作 $deg A$

对于多项式 $A(x), B(x)$ ，存在唯一的 $Q(x), R(x)$ 满足 $A(x) = Q(x)B(x) + R(x)$ ，其中 $deg R < deg B$ ，我们称 $Q(x)$ 为 $B(x)$ 除 $A(x)$ 的商， $R(x)$ 为 $B(x)$ 除 $A(x)$ 的余数，可以记作

$A(x) \equiv R(x) \pmod {B(x)}$

(more…)

从多项式乘法到快速傅里叶变换

April 25, 2015January 21, 2018 miskcoo 31 Comments

概述

计算多项式的乘法，或者计算两个大整数的乘法是在计算机中很常见的运算，如果用普通的方法进行，复杂度将会是 $\mathcal O(n^2)$ 的，还有一种分治乘法，可以做到 $\mathcal O(n^{\log_23})$ 时间计算（可以看这里）。下面从计算多项式的乘法出发，介绍快速傅里叶变换（Fast Fourier Transform, FFT）如何在 $\mathcal O(n\log n)$ 的时间内计算出两个多项式的乘积。另外，存在只需要两次快速傅立叶变换就可以计算大整数乘法的方法，具体见实序列离散傅里叶变换的快速算法

准备知识

这里介绍一些后面可能会用到的知识（主要是关于多项式、卷积以及复数的），如果你已经知道觉得它太水了或者想用到的时候再看就跳过吧

多项式

简单来说，形如 $a_0+a_1X+a_2X^2+\cdots+a_nX^n$ 的代数表达式叫做多项式，可以记作 $P(X)=a_0+a_1X+a_2X^2+\cdots+a_nX^n$ ， $a_0, a_1, \cdots, a_n$ 叫做多项式的系数， $X$ 是一个不定元，不表示任何值，不定元在多项式中最大项的次数称作多项式的次数

多项式的系数表示法

像刚刚我们提到的那些多项式，都是以系数形式表示的，也就是将 $n$ 次多项式 $A(x)$ 的系数 $a_0, a_1, \cdots, a_n$ 看作 $n+1$ 维向量 $\vec a=(a_0,a_1,\cdots,a_n)$ ，其系数表示（coefficient representation）就是向量 $\vec a$

多项式的点值表示法

如果选取 $n+1$ 个不同的数 $x_0, x_1, \cdots, x_n$ 对多项式进行求值，得到 $A(x_0), A(x_1), \cdots, A(x_n)$ ，那么就称 $\{\left(x_i, A(x_i)\right) : 0 \leq i \leq n, i \in \mathbb Z\}$ 为多项式 $A(x)$ 的点值表示（point-value representation）

多项式 $A(x)$ 的点值表示不止一种，你只要选取不同的数就可以得到不同的点值表示，但是任何一种点值表示都能唯一确定一个多项式，为了从点值表示转换成系数表示，可以直接通过插值的方法

复数

后面提到的 $i$ ，除非作为 $\sum$ 求和的变量，其余的都表示虚数单位 $\sqrt{-1}$

单位根

$n$ 次单位根是指能够满足方程 $z^n=1$ 的复数，这些复数一共有 $n$ 个它们都分布在复平面的单位圆上，并且构成一个正 $n$ 边形，它们把单位圆等分成 $n$ 个部分

根据复数乘法相当于模长相乘，幅角相加就可以知道， $n$ 次单位根的模长一定是 $1$ ，幅角的 $n$ 倍是 $0$

这样， $n$ 次单位根也就是

$e^{\frac{2\pi ki}{n}}, k = 0, 1, 2, \cdots, n - 1$

再根据欧拉公式

$e^{\theta i}=\cos\theta + i\sin\theta$

就可以知道 $n$ 次单位根的算术表示

如果记 $\omega_n=e^{\frac{2\pi i}{n}}$ ，那么 $n$ 次单位根就是 $\omega_n^0, \omega_n^1, \cdots, \omega_n^{n-1}$

多项式的乘法

给定两个多项式 $A(x), B(x)$

$A(x) = \sum_{i=0}^na_ix^i = a_nx^n+a_{n-1}x^{n-1}+\cdots+a_1x+a_0 \\ B(x) = \sum_{i=0}^nb_ix^i = b_nx^n+b_{n-1}x^{n-1}+\cdots+b_1x+b_0$

将这两个多项式相乘得到 $C(x) = \sum_{i=0}^{2n}c_ix^i$ ，在这里

$c_i=\sum_{j+k=i,0\leq j,k\leq n}a_jb_kx^i$

如果一个个去算 $c_i$ 的话，要花费 $\mathcal O(n^2)$ 的时间才可以完成，但是，这是在系数表示下计算的，如果转换成点值表示，知道了 $A(x), B(x)$ 的点值表示后，由于只有 $n+1$ 个点，就可以直接将其相乘，在 $\mathcal O(n)$ 的时间内得到 $C(x)$ 的点值表示