FFT与多项式乘法

正大光明松鼠

已于 2024-09-26 18:48:28 修改

阅读量1.1k

点赞数 27

文章标签：密码学区块链

于 2024-09-04 17:33:47 首次发布

本文链接：https://blog.csdn.net/sunguizhu/article/details/141860677

版权

本文参考文献：

Polynomial Multiplication (Vector Convolution) and FFT

多项式 I：拉格朗日插值与快速傅里叶变换 - 洛谷专栏

前言

通常两个多项式相乘所需的时间复杂度是 Θ(n^2)，使用FFT算法可以将多项式乘法的时间复杂度降为 $\Theta (nlog_2n)$ 。

1.多项式表示方法

这节介绍两种表示多项式的方法：系数表达和点值表达。

1.3.1. 系数表达

对一个次数界为 n 的多项式 $f(x) = a_{0} + a_{1}* x + ...+ a_{n-1}* x^{n-1}$ 而言，其系数表达是一个由系数组成的向量 $a = (a_{0}, a_{1}, ..., a_{n-1})$ 。

采用系数表达对于多项式的某些运算是很方便的。例如，对多项式 f(x) 在给定点 $x_{0}$ 的求值运算就是计算 $f(x_{0})$ 的值。使用霍纳法则，我们可以在 Θ(n) (n-1次乘法和n-1次加法)时间复杂度内完成求值运算：

$f(x_0) = a_0 + x_0(a_1+x_0(a_2+...+x_0(a_{n-2} + x_0(a_{n-1})...))$

1.3.1.1多项式的运算

假设有多项式f(x)和g(x)：

$f(x)$ 的系数是 $a = (a_{0}, a_{1}, ..., a_{n-1})$ ,

$g(x)$ 的系数是 $b = (b_{0}, b_{1}, ..., b_{n-1})$ ,

1.加法运算：

在系数表示形式下，求两个多项式的加法是很容易的。

令 $h(x) = f(x) + g(x)$ ，假设 $h(x)$ 的系数是 $c = (c_{0}, c_{1}, ..., c_{n-1})$ ，显然有 $c_i = a_i + b_i$

因此计算 $h(x)$ 的计算复杂度是 Θ(n) 。

2.乘法运算

令 $h(x) = f(x) * g(x)$ ，则 $h(x) = (\sum_{i=0}^{n-1} a_{i}x_{i})*(\sum_{j=0}^{n-1} b_{j}x_{j})=\sum_{i=0}^{2n-2}(\sum_{j=0}^{i}a_jb_{i-j})*x_{i}$ ,

因此计算复杂度是 Θ(n^2) (需要进行 $(n-1)^2$ 次乘法)

本文的主要目标就是 讨论如何降低多项式乘法的时间复杂度。

1.3.2. 点值表达

一个次数界为 n 的多项式 f(x) 的点值表达就是一个由 n 个不同的点值对组成的集合：

${(x_0,y_0), (x_1,y_1),...,(x_{n-1},y_{n-1})}$ , 其中 $y_k = f(x_k)$ ，且 $x_k$ 各不相同。

1.3.2.1 点值表达式的运算

假设已知f(x)的n-1个点： ${(x_0,y_0), (x_1,y_1),...,(x_{n-1},y_{n-1})}$

假设已知g(x)的n-1个点： ${(x_0,y_{0}^{'}), (x_1,y_{1}^{'}),...,(x_{n-1},y_{n-1}^{'})}$

由霍纳法则知，求多项式的一个点需要的时间是n,则求n个点需要的时间是n*n，即时间复杂度是 $\Theta (n^2)$

1.加法运算：

令 $h(x) = f(x) + g(x)$ ，则 $h(x)$ 的n个点为：

${(x_0,y_0 + y_{0}^{'}), (x_1,y_1 + y_{1}^{'}),...,(x_{n-1},y_{n-1} + y_{n-1}^{'})}$

计算 $h(x)$ 的计算复杂度是 Θ(n) 。

2.乘法运算

令 $h(x) = f(x) * g(x)$ ，

由于h(x)的次数是n-2,因此f(x)和g(x)需要各选取2n-1个点,

假设已知f(x)的2n-1个点： ${(x_0,y_0), (x_1,y_1),...,(x_{2n-1},y_{2n-1})}$

假设已知g(x)的2n-1个点： ${(x_0,y_{0}^{'}), (x_1,y_{1}^{'}),...,(x_{2n-1},y_{2n-1}^{'})}$

则h(x)的点值表达为： ${(x_0,y_{0}y_{0}^{'}), (x_1,y_{1}y_{1}^{'}),...,(x_{2n-1},y_{2n-1}y_{2n-1}^{'})}$

计算h(x)计算复杂度是 Θ(n) 。

1.3.3. 系数表达和点值表达的相互转换

从系数表示法转为点值表示法称为求值（Evaluation）。
从点值表示法转为系数表示法称为插值（Interpolation）。

1.3.3.1 求值法

左侧矩阵就是著名的 范德蒙德矩阵。

当 m=n时为范德蒙德方阵， $x_i$ 互不相同时其逆存在，帮助我们快速从点值表示法转回系数表示法。m>n时任取 xi 互不相同的 n+1 行可以求逆。m<n 时无法还原系数。这体现出 n+1个点值唯一确定最高次数不超过 n 的多项式。

朴素计算求值的复杂度为 O(nm)，因为带入求值一次的复杂度为 O(n)。快速傅里叶变换即在离散傅里叶变换基础上通过选取合适的 xi，使得可以快速求值。

1.3.3.2拉格朗日插值法

从一个多项式的“点值表达”确定其“系数表达”形式称为“插值”。目前最常使用的插值方法即为拉格朗日插值法。

假设已知f(x)的n-1个点： ${(x_0,y_0), (x_1,y_1),...,(x_{n-1},y_{n-1})}$ 。则根据拉格朗日插值定理知，f(x)的多项式如下：

$f(x) = \sum_{i=0}^{n-1}y_i\prod_{j=0,j\neq i}^{j=n-1}\frac{x-x_j}{x_i-x_j}$

扩展：

我们称 $L_x(x) = \prod_{j=0,j\neq i}^{j=n-1}\frac{x-x_j}{x_i-x_j}$ 为拉格朗日基，其特点是:

当且仅当 $x = x_i$ 时，值为1， $x$ 为其他值时，值为0。

我们可以用这一性质对多项式进行约束。

假设某个多项式a(x)在 $x_0$ 处的值是 $y_0$ ，则要求a(x)满足下面等式：

$L_0(x)*(a(x) - h) = 0$ ，对于任意 $x\in (x_0,x_1,...x_{n-1})$ 成立

1.3.4 总结

根据上面分析，可以发现两种形式的多项式乘法，算法复杂度均为 $\Theta(n^2)$ 。

可以看出点值表示法的大部分计算时间在于求值和插值方法。如果这两部分可以提速，那么点值法的计算复杂度变可以降低。

很幸运，我们通过在特殊点对多项式进行取值，可以将求值和插值部分,分别使用FFT和IFFT运算，从而将计算复杂度变降为 $\Theta (nlog_2n)$ 。

2.单位根

由于FFT会用到单位根，因此在介绍FFT之前，先介绍一下单位根。

如果 $g^n = 1$ ，则g成为n次单位根。本定义通常用于复数域和有限域。

下面的介绍均为复数域范畴。

欧拉公式： $e^{ix} = cos x + i sin x$ .

(cos为点在x轴上的坐标，sin为点在y轴上的坐标)

即单位圆上从 (1,0)开始旋转 x度得到的复数。

由于 $e^{i*2\pi } = cos2\pi +isin2\pi =1 +i * 0 = 1$ .

故有 $(e^{i*\frac{2\pi }{n} } )^n = e^{i*2\pi }= cos2\pi +isin2\pi =1 +i * 0 = 1$

故 $e^{i*\frac{2\pi }{n} }$ 是n次单位根，可记为 $w_n$ 。另外 $w_{n}^{k} = e^{i*k\frac{2\pi }{n} }, k\epsilon (0,n-1)$ 均为n次单位根。

由于 $cos \pi = -1, cos 2\pi = 1, sin \pi = sin 2\pi = 0$ ,

故 $w^{k+\frac{n}{2}}_{n} = e^{i*(k+\frac{n}{2})\frac{2\pi }{n} } = e^{i*(k\frac{2\pi }{n} +\pi )} = e^{i*k\frac{2\pi }{n}} *e^{i\pi } = e^{i*k\frac{2\pi }{n}} * (-1) =-w_{n}^{k}$

故单位根有如下性质：

3.FFT

3.1 DFT

在介绍 FFT 之前，我们先给出离散傅里叶变换（Discrete Fourier Transform，DFT）的概念。

DFT 在工程中是将离散信号从时域转为频域的过程。碰巧的是，其表达式刚好可以用来对多项式进行多点求值，只不过这些点值是固定的 单位根 处的点值，但对于求值做多项式乘法已经足够了。

DFT即为对f(x)在n个n次单位根处求值。如下：

设 $f(x) = a_{0} + a_{1}* x + ...+ a_{n-1}* x^{n-1}$ ，

则 $y_k = f(w_{N}^{-k}) = \sum_{i=0}^{n-1} a_i(w_{N}^{-k})^n$

3.2 FFT算法详解

3.2.1 简化情况

首先我们得搞清楚，DFT 是一个变换，而 FFT 是用于实现 DFT 的算法。

FFT的巧妙之处在于通过将函数拆分为奇函数和偶函数，将函数的次数降低了，如此可以迭代 $log_{2}^{n}$ 次，最终将算法复杂度降低到 $\Theta(nlog_{2}^{n})$

下面将具体讲解这个过程。

任意多项式f(x)可以拆分为偶函数 $f_e(x)$ 和奇函数 $f_o(x)$ 之和，则：

$f_{e}^{'}(u)$ 和 $f_{o}^{'}(u)$ 的系数均n/2个,次数也为n/2，分别计算 $f_{e}^{'}(u)$ 和 $f_{o}^{'}(u)$ 在n/2个点处的值，就可以得到f(x)在n个点上的取值（因为可以同时计算出f(x)和f(-x)两个点）。

因此对于这两个多项式中的任意一个求一个点的复杂度为 $\Theta(\frac{n}{2})$ ,求n/2个点的复杂度是 $\Theta(\frac{n}{2} * \frac{n}{2}) = \Theta((\frac{n}{2})^2)$ 。则对2个多项式均求n/2个点的复杂度是 $2 * \Theta((\frac{n}{2})^2)$ 。

3.2.2 引入单位根

上面的式子拆分奇偶可以一直迭代的关键是，每次迭代都要找到互为相反数的两个x坐标。但是每次迭代，新的变量u都是原变量x的平方，在实数域中，平方是没有复数的。

由于单位根始终有如下性质： $w^{k+\frac{n}{2}}_{n} = e^{i*(k+\frac{n}{2})\frac{2\pi }{n} } = e^{i*(k\frac{2\pi }{n} +\pi )} = e^{i*k\frac{2\pi }{n}} *e^{i\pi } = e^{i*k\frac{2\pi }{n}} * (-1) =-w_{n}^{k}$ ，

故始终有互为相反数的一组数据存在，所以此时可以引入单位复根。

因此，FFT的范德蒙方阵如下：