FFT实现大整数乘法供参习.docVIP

目录 FFT实现高精度乘法 3 绪论 3 DFT（离散傅里叶变换） 5 Cooley-Tukey算法实现FFT（快速傅里叶变换） 5 位元反转（Bit reversal） 8 C++算法实现 9 参考资料： 12 FFT实现高精度乘法 绪论 对于这个课题，我们先从多项式开始一步一步进行分析。一个度数为n的多项式A(x)可以表示为： 一个多项式可以有不同的表示方式，一种是我们常见的系数表示，另一种是点值表示。例如上面的多项式就是以系数表示形式定义的。它展开之后可以写成下面这样： 这种方式的优点是很直观，也是我们最常用的表示多项式的方法，但是采用这种方式表示的多项式在进行大数值的乘法运算的时候，如果采用直接乘法方式，假设两个多项式的度数都为n, 那么乘法计算的时间复杂度为O(n2). 如果采用另一种不那么直观的点值表示法，对于多项式的度数比较大的情况，采用点值表示法进行乘法运算可以获得较大的性能提升。对于上面提到的多项式A(x), 它的点值表示法大概像下面这样子： 其中x0~xn-1是A(x)上面取的不同的点，而且 可以看出，多项式从系数表示形式直接转化为点值表示形式需要选取n（多项式的度数）个不同的值，进行n次的带入求值运算。那么转化为点值表示形式有什么好处呢？好处在于两个满足一定要求的点值表达式进行乘法运算的时间复杂度为O(n). 相对于系数表示法直接进行乘法运算，点值表达式在计算乘法的性能是很理想的。 现在我们来看看点值表达式是如何快速地进行多项式乘法的。假设有两个多项式A(x), B(x)，它们的乘积为一个新的多项式C(x)。在计算乘法的时候，多项式的项数可能会增加，因此乘法运算得到的新的多项式C(x)的度数是degree(C) = degree(A) + degree(B). 因为这个原因，在计算乘法的时候，需要对A(x)和B(x)进行点值表达式的扩张。也就是说多项式A(x), B(x)的点值表示形式要分别扩张为： 和 还有一点要注意的是上面两个式子中的x0~x2n-1是要一一对应的，否则不能直接进行乘法运算。这样C(x)就可以转化为A(x), B(x)的点值表达式的对应项相乘的结果： 这样，我们就能以O(n)的时间复杂度（不考虑点值表达式的扩张运算，只需要计算2n次乘法）进行多项式A(x)和B(x)之间的乘法运算。 那么接下来，我们要讨论的是如何进行多项式的点值表示形式和系数表示形式之间的转换，这里涉及到一个问题，就是多项式的插值转换的唯一性。也就是说你要确定从点值形式转化为系数形式是唯一的，这个问题已经被数学证明，插值转换是唯一的，这个是利用线性代数的范德蒙德矩阵（下图左方的矩阵）的可逆性来证明的，具体的证明就不贴上来了。而系数形式转化为点值形式可以有很多种不同的方案，因此不存在唯一性。 对于系数表达式转化为点值表达式，[1]称它为求值(evaluate)，我们可以利用秦九韶算法对给定的n个点进行代入多项式A(x)运算，这种方法的时间复杂度为O(n2). 对于点值表达式转化为系数表达式，[1]称它为插值(interpolate)，我们可以直接利用上图给出的线性代数方程，求出范德蒙德矩阵的逆矩阵，然后与向量上图右方的y向量相乘，计算出系数向量a, 用这种方法计算的时间复杂度从上图就能看出来，为O(n3). 也可以利用拉格朗日公式进行插值运算。 按拉格朗日公式计算系数表达式的时间复杂度为O(n2). 采用以上提到的方法进行系数表示形式和点值表示形式之间的转化的效率并不理想，因此我们要考虑如何快速进行两者之间的转化，这就是FFT需要解决的问题。 以下这张图就是解决这个问题的总体概览，我们不妨一窥全豹： DFT（离散傅里叶变换） 上图中的求值(evaluation)和插值(interpolation)是分别通过DFT和反向DFT实现的。DFT是傅里叶分析的一种，它将一组有限的，均匀分布的函数样本值转化为一组复变正弦函数的系数。因为我看的资料大部分都是英文的，有些专业的词汇好像没有对应的翻译，为了严谨，这部分我就简单地说说DFT是怎么做的，尽力避免陈述其它的数学知识。 假定我们要进行对度数为n的多项式A(x)进行DFT，我们要将下列的值 代入A(x)分别求值，得到一个y向量，这个向量就是原先多项式的系数向量的DFT，可以写为。 上式出现的叫做principal nth root of unity, 它的值为 这是最原始的计算DFT的方法，按这种方法计算的话，时间复杂度是O(n2). 接下来我们讨论如何进行反向DFT，其实反向DFT是很简单的，根据我们上面提到的插值转换的唯一性理论(Uniqueness of an interpolating polynomial)，我们