大数乘法（快速傅立叶变换）上

在谈论大数乘法前，先来看看多项式乘法，假设

A = 7x^3 + 5x^2 + 3x + 4

B = 4x^2 + 6

C = A * B

那传统的做法就是这样：

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 7x^3 + 5x^2 ?+ ?3x ?+ 4

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?4x^2 ? ? ? ? ? ?+ 6

? ? ? ? ? ? ? ? ? ? ? ?---------------------------------------

? ? ? ? ? ? ? ? ? ? ? ? ? ? +42x^3 + 30x^2 + 18x + 24

? ? ?28x^5 + 20x^4 +12x^3 + 16x^2

-----------------------------------------------------------

? ? 28x^5 + 20x^4 + 54x^3 + 46x^2 + 18x + 24

若A 和 B 都有 n 个系数，很明显这样做的时间复杂度为O(n^2)，当n很大时效果相当不理想

? ? ? ? 像?B = 4x^2 + 6 这样的式子叫系数表示法；如果我们把 B 一般化，即 B = b0 + b1x + b2x^2；则我们至少需要 3 对 (x,y) 值才能求出 b0,b1,b2。这里 b0 = 6,b1 = 0,b2 = 4；例如 B 在x0 = 1,x1 = 2,x2 = 3 的情况下的解分别为 y0 = 10,y1 = 22,y2 = 42。那我们也可以用 {(1,10)，(2,22)，(3,42)}，即 {(x0,y0)，(x1,y1)，(x2,y2)} 来表示唯一的多项式B，即假如我们在不知道b0,b2的情况下可以用{(1,42)} 来唯一表示多项式B，即唯一确定b0,b2。当然 x0,x1,x2 也可以取其他值，但对于这个例子一定要取够3个，因为有3个未知数 b0,b2；对应的 y0,y1,y2 也会不同，这就叫点值表示法。

? ? ? ? 为什么要这样做？还是上面那个例子，我们用点值表示法算一遍C = A * B，（A 和 B 的xi要对应相同）还是假设取x0 = 1,x2 = 3，我们发现 A 有4个未知数，所以要多一个x3，假设x3 = 4；为了方便计算，我们把 B 的 (x3,y3) 也算出来。

A = {(1,19)，(2,86)，(3,247)，(4,544)}

B = {(1,42)，(4,70)}

xi 不变，yi对应相乘得到 Ci = (xi,yai * ybi)

C =?{(1,190)，(2,1892)，(3,10374)，(4,38080)}

到这里我们发现对于一般情况的等长度的两个多项式相乘，结果的长度会接近加倍，比如

A = a0 + a1x + a2x^2 + a3x^3

B = b0 + b1x + b2x^2 + b3x^3

C = A * B = c0 + c1x + c2x^2 + c3x^3 +?c4x^4 + c5x^5 + c6x^6

? ? ? ? ?因此有n个系数的两个多项式相乘（比如A * B）时仅用其n个点值相乘是不能正确表示结果的（比如C），所以我们要使用足够多的点值，比如这里 A 和 B 都应该至少使用7个点值。如下（数有点大，偷下懒。。。）

A = {(x0,y2)，(x3,y3)，(x4,y4)，(x5,y5)，(x6,y6)}

B = {(x0,y0')，(x1,y1')，(x2,y2')，(x3,y3')，(x4,y4')，(x5,y5')，(x6,y6')}

C = {(x0,y0*y0')，(x1,y1*y1')，(x2,y2*y2')，(x3,y3*y3')，(x4,y4*y4')，(x5,y5*y5')，(x6,y6*y6')}

? ? ? ? 这样C的点值才能表示唯一的6次多项式（7个未知系数有6次幂），最后，用某方法（重点下面讲）把C的点值表示法转换成系数表示法就达到了目的，这个过程叫插值

? ? ? ? 所以，先理下思路

1）分别计算A、B足够多的点值（也叫求值），使其能唯一表示C（如A、B的系数都有n个，则通常扩展到2n个点值）

2）点乘：即A、B点值对应相乘得到C

3）插值：把C的点值表示法转换为系数表示法

? ? ? ? 其中点乘的复杂度显然为O(n)，我们记N = 2n（记住了下面要用的~），所以我们如能找到一种方法能够快速对A、B进行求值，插值，那目的就达到了，事实上也的确有这样的方法，这也是多项式乘法的重点，不然都没法往下讲。。。

?----------------------------------------------------------------------------------------------------------------------------------------------------

? ? ? ? 先看求值，y = a0 + a1x + a2x^2 + ... + anx^n ?给出一个 x 要求 y 如果按正常顺序求需要做 (1 + n) * n / 2 次乘法和 n 次加法；这种方法肯定是不可取的，因为算一个值都要 O(n^2)，那算 N 个值就是 O(n^3) 了，比系数表示法直接乘还慢~

? ? ? ? 再来看另一种方法，国产方法，没错就是它，秦九韶算法（在国外也叫霍纳(Horner)法则）；为了理解方便，我们举个短的例子 y = a0 + a1x + a2x^2 + a3x^3；那它就可以化为 y = a0 + x(a1 + x(a2 + a3x))，这样对某个x求y就需要做 n 次乘法和 n 次加法，算 N 个值的复杂度为 O(n^2)。比上面的好了点，不过跟系数表示法相比并没有什么优势~

? ? ? ? 所以我们需要第三种方法，这种方法的特点其实是在 xi 的选取上很特殊；它选的 xi 不是从 x0 = 1,...,xn = n 或者其他一些相异的整数或实数；它选的是复数！当然也不是乱选的，因此我们要先在这里停一下，吸口气，去打下复数的基础。。。

---------------------------------------------------------------------------------------------------------------------------------------------------------

? ? ? ? ?

? ? ? ? ?z = x + yi 这是常见的复数表示形式大家都知道，|z| = sqrt(x^2 + y^2) 表示复数的模

? ? ? ? 用建立了笛卡尔直角坐标系的平面表示复数的平面叫做复平面

? ? ? ? 下面必须盗图了（来自《复变函数与积分变换》）不难的，请坚持看完！

? ? ? ? 复数的指数式也称为欧拉公式，它更像是一种定义。欧拉公式是非常美妙的式子，为什么要引入它呢？主要是为了简化计算！因为 z = x + yi 可以表示成 z = r*(cosθ + i*sinθ)，r = |z|。看个例子：

?

? ? ? ? 所以用指数式表示复数是很简洁的，我们最后的求值与插值的方法也是用复数的指数式表示。

? ? ? ? 这里有一点要提一下，有人说?θ = 0 时，e^(θi) ?=?e^(0*i) = 1；结果没有错，但不是这样算的。e^(θi) 不能理解为?e^(θ*i)；那个不是乘号，只是一种记号。正确解法应该把它展开算?e^(θi) ?= cosθ + isinθ =?cos0?+ isin0 = 1

? ? ? ? 设 z 为复数，n 为正整数，若存在复数 w 满足方程 w^n = z；

? ? ? ? ?重要的事情说三遍：

? ? ? ? ?在几何上 z^(1/n) 的 n 个值就是以原点为中心，r^(1/n) 为半径的圆的内接正 n 边形的 n 个顶点。

? ? ? ? ?在几何上 z^(1/n) 的 n 个值就是以原点为中心，r^(1/n) 为半径的圆的内接正 n 边形的 n 个顶点。

? ? ? ? ?在几何上 z^(1/n) 的 n 个值就是以原点为中心，r^(1/n) 为半径的圆的内接正 n 边形的 n 个顶点。

? ? ? ? ?还有上面的公式 z = e^(iθ) != 0，wk = (r^(1/n)) * e^(i(θ+2*k*pi) / n)。θ 为 z 在复平面中向量表示的幅角。

? ? ? ??

? ? ? ? ?做个题目理解下：

? ? ? ? 这个就不画了，想加深理解可以自己画一下~

--------------------------------------------------------------------------------------------------------------------------------------------------------

? ? ? ? 接下来才是高能部分：

? ? ? ? 我们前面讲了，在对 A = a0 + a1x + a2x^2 + ... + a(N-1)x^(N-1) （我们假设A的系数有 N（前面讲过N = 2n） 个，即a0 ~ a(N-1)）求值时我们选的 N 个 xi (0 <= i < N) 用的是复数，那现在告诉你这 N 个 xi 就是 w^N = 1 的 N 个解，w0 ~ wN 叫N 次单位复根。即 x0 ~ xN 分别等于 w0 ~ wN；

? ? ? ? 这里的 N 也就是一种记号，如果看不顺，理解成 n 也是没问题的~

? ? ? ? 再上一遍公式

? ? ? ? 对于 N 次单位复根的求解，θ 值为 0 ，r^(1/N) = 1。其中解 w1 = e^(i(2*pi/N)) 叫主单位复根；其实也就是把一个圆分成 N 份，w1 表示每份的大小，只不过是在复平面上表示。

? ? ? ? 下面我们要换下记号了

? ? ? ? ?记 N 为 n

? ? ? ? ?则 A = a0 + a1x + a2x^2 + ... + a(n-1)x^(n-1)

? ? ? ? ?也可记为?

? ? ? ? ?原来的各个单位复根 w0,w1,w(n-1) 记为?

，下标的 n 就是对单位复数 1 开 n 次方的意思，很简单大家都能明白。这样是因为下面的计算会改变下标。

? ? ? ? ?那用?

?代入 xk 来求解 A(xk) 就可以表示为?

? ?(0 <= k <= n-1)

? ? ? ? ?如果把 A 的系数用向量表示 a = {a0,a1,a2,a(n-1)}，把求出的所有 y 值也用向量表示 y = {y0,y2,y(n-1)}，那我们可以称向量 y 为 向量 a 的离散傅立叶变换（Discrete Fourier Transform,DFT）。我们的目的之一就是把向量 a 的离散傅立叶变换（即向量 y ）快速地算出来，我们使用的是一种叫快速傅立叶变换（FFT）的算法。这里我们可以猜到，FFT 算法起码要比系数表示法直接相乘的 O(n^2) 的复杂度更低才行，幸运的是，FFT 的时间复杂度是 O(n*logn)；下面讲的是 FFT 是怎样快速地计算向量 y 的。

? ? ? ? ?在讲FFT之前需要讲清楚一件事，就是上面所讲的 n 到底是多大？

? ? ? ? ?举个例子吧，比如：

? ? ? ? ?A = 3 + 2X^2 + 4X^3 ?+ 6x^4

? ? ? ? ?B = 2x + X^2 + 5X^3

? ? ? ? ?以 A 举例（B类似），它可以写成标准式：A = a0 + a1x + a2x^2 + a3x^3 + a4x^4；这里有 5 个系数，根据前面讲的为了能正确表示 A 与 B相乘的结果 C，我们把 A 放大到 A = a0 + a1x + ... + a9x^9，使它变成10个系数（当然 a5 ~ a9 的值都是 0，B 类似），同时用 A 和 B 的 10 个点值计算C。不过对于 FFT 算法，这样扩大是不行的；根据算法要求（下面会讲清楚，现在不用管），我们要确保 n 为 2 的幂，所以 A 有 5 个系数，不是 2 的幂，我们应该把它先扩大到比5大且离 5 最近的那个 2 的幂，即 2^3 = 8。然后再把 8 扩大两倍变成 16；这才是 A B 两个多项式相乘的最终的 n。当然原来系数个数就是 2 的幂的话就只需要把它扩大两倍就行了。当系数个数为 513 的话，没错，要先把它扩大成1024个系数，然后再放大两倍变成 2048个系数；就是这么狠！想想这样有点耗空间~

? ? ? ? ?

? ? ? ? ?所以最终的 A ?= a0 + a1x + a2x^2 + a3x^3 + ... + a(n-1)x^(n-1) 中的系数个数，即 n 必须是 2 的幂，我们总能够通过添加 ai = 0 满足需要。

? ? ? ? FFT 算法运用了分治策略，这下有点理解为什么 n 要是 2 的幂了吧~，什么是分治策略就不展开讲了，不懂的就搞清楚了再往下看。

? ? ? ??

? ? ? ? 次数界就是多项式最高次幂是多少次方。

? ? ? ? 莫慌，我们看看 A(k) = a0 + a1k + a2k^2 + a3k^3 + ... + a(n-1)k^(n-1)?怎样算

? ? ? ? 把 k^2 代入?

? 和

?有

? ? ? ? ? ?

? ? ? ? 不难看出 A(k) =?

?+ k *?

? ? ? ? 这样，求次数界为 n 的多项式?A(x) 在?

处的值就转化为：

? ? ? ? 求次数界为 n / 2 的多项式?

?和?

?在?

?处的值。

? ? ? ? 形象地举个 n = 4 的例子就是：

? ? ? ?

? ? ? ? 在这里我们引入定理：

? ? ? ? 相消引理：对任何整数 n >= 0，k >= 0，d > 0；?

? ? ? ? ? ? ? ? ? ? ? ? ? ?也不难证明：

? ? ? ? ? ? ? ? ? ? ? ? ? ?可以得出推论：对于任意偶数 n > 0， 有 ?

（举个例子上下标不断除以2嘛）

? ? ? ? 结合上面的 n 个点??

，其实只有 n/2 个不同的值。

? ? ??? 上面 n = 4 的例子中

? ? ?

? ? ? ? 可以知道?

，因此?

? ? ? ? 所以，??

前面 n/2 项跟后面 n/2 项是重复出现的；如上面?

?和?

?重复出现了2次。

? ? ? ? 因此，之前的式子可以改为：

? ? ?

? ? ? ?这样，求次数界为 n 的多项式?A(x) 在?

处的值就完全转化为：

? ? ? ?求次数界为 n / 2 的两个多项式?

?和?

?在?

?处的值；规模神奇地缩小了一半，之后我们递归处理

和?

就可以把向量 y 整个儿求出来！伪代码应运而生：

? ? ? ? 如果伪代码还需要解释的话就说明之前的知识没有理解透彻，请往回多看几次（主要是高能部分）！

? ? ? ? 这就是多项式系数表示法转换为点值表示法（也称为求值过程）的 O(n*logn) 复杂度所用到的FFT算法的整个过程。

? ? ? ? 所以我们可以将多项式 A 和 B 按照上面的方法分别转化为点值表示，然后 y 值对应相乘就得到了多项式 C 的点值表示，接下来把 C 的点值表示转换为系数表示就完成了多项式乘法的运算。所以，下面讲介绍插值过程的算法。

? ? ? ? 插值过程的算法：?大数乘法（快速傅立叶变换）下