FFT实现大整数乘法

目录FFT实现高精度乘法3绪论3DFT（离散傅里叶变换）5Cooley-Tukey算法实现FFT（快速傅里叶变换）5位元反转（Bit reversal）8C+算法实现9参考资料：12FFT实现高精度乘法绪论对于这个课题，我们先从多项式开始一步一步进行分析。一个度数为n的多项式A(x)可以表示为：一个多项式可以有不同的表示方式，一种是我们常见的系数表示，另一种是点值表示。例如上面的多项式就是以系数表示形式定义的。它展开之后可以写成下面这样：这种方式的优点是很直观，也是我们最常用的表示多项式的方法，但是采用这种方式表示的多项式在进行大数值的乘法运算的时候，如果采用直接乘法方式，假设两个多项式的度数都为n, 那么乘法计算的时间复杂度为O(n2).如果采用另一种不那么直观的点值表示法，对于多项式的度数比较大的情况，采用点值表示法进行乘法运算可以获得较大的性能提升。对于上面提到的多项式A(x), 它的点值表示法大概像下面这样子：其中x0xn-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)的点值表示形式要分别扩张为：和还有一点要注意的是上面两个式子中的x0x2n-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)，我们进行DFT的运算实质是下面这样：这样，我们只要求出中间的范德蒙德矩阵的逆矩阵，然后将与向量相乘，就能得出系数向量，这就是反向DFT，计算的时间复杂度为O(n3)。能这样计算的前提是上图的范德蒙德矩阵是可逆的，当然它的可逆性已经被证明了，我这里就不写了。从上面的讨论可以看到无论是DFT还是反向DFT，直接计算的时间复杂度都不是很理想，接下来我们就来讨论一下FFT是如何降低计算DFT的时间复杂度的。Cooley-Tukey算法实现FFT（快速傅里叶变换）上文提到了DFT是什么以及如何进行DFT和反向DFT运算，现在我们来讨论FFT是如何快速进行DFT的。我们在上文提到过一个数值，这个数值有个有意思的特性，根据Halving lemma, 下式是成立的：其中k为非负整数。这个等式是FFT得以实现的关键。值得一提的是FFT的实现的算法有多种，我们这里讨论的FFT算法是基于1的。这本书上提到的FFT算法叫做Cooley-Tukey算法，它是目前为止应用得最广泛的FFT算法。从现在开始，如果没有别的说明，我们就以FFT算法来称呼Cooley-Tukey算法。假设我们有一个以系数形式表示的多项式A(x)，它的度数为n，n为2的幂。那么它可以表示为现在我们将它分成更小的两个多项式和，它们分别满足如下关系：那么拆分出来的两个多项式和原来的多项式有如下关系：这样的处理就是分治的基本步骤，将代入进行DFT的问题就可以转化为：1、将代入度数为n / 2的和分别求值；2、将计算出来的结果根据上文提到的递归关系合并成最终的结果。根据Halving lemma，并不是n个各不相同的值，它有一半是重复的。这样我们就能将问题的规模缩小成原来的二分之一，这是能够利用递归提升效率的关键。现在我们就能够根据上面的分析写出FFT的递归算法了，但递归算法的效率不高，接下来我们要讨论FFT算法的非递归实现。我们来分析一下递归形式的FFT的步骤是怎么样的，假设我们要对元素个数为8的向量进行FFT运算，下面是自顶向下的递归的步骤：运算的逻辑步骤是一颗树，而且我们猜想叶子结点的元素排列序号与根节点的元素排列序号应该有一定的联系。我们只要找出两者之间的联系，就可以通过它们之间的关系直接计算出递归最深层的元素排列序列，然后直接从最底层开始进行运算，用计算出的DFT替换掉上一层的节点，重复这个步骤直到根节点为止。这样我们就能将自顶向下的递归形式的FFT算法改进成自底向上的迭代形式的算法了。为了分析它们之间的关系，我们打算从元素的下标着手，原始序列是0, 1, 2, 3, 4, 5, 6, 7, 它们的二进制表示分别为000, 001, 010, 011, 100, 101, 110, 111. 最底层的元素排列序列的下标是0, 4, 2, 6, 1, 5, 3, 7. 它们的二进制表示分别为000, 100, 010, 110, 101, 011, 111. 以下的表格可以清楚地显示对应位置的变化:原始序列000001010011100101110111最终序列000100010110001101011111从表中可以看出原始的序列和最终的序列的对应位置的下标是互为相反位。也就是说原始序列的对应位置的下标进行位元反转运算就能得到最终序列。这样我们就能够利用这个关系着手实现FFT算法的迭代形式了。算法实现的伪代码如下：其中BIT-REVERSE-COPY将我们要进行FFT的a向量按照位元反转位置复制到另一个向量A，复制完成之后的A就是我们要进行递归运算的最终序列。它的伪代码如下：其中rev实现位元反转，关于位元反转的具体问题我们在后面还会提到，现在我们来分析一下FFT非递归算法解决问题的步骤。在ITERATIVE-FFT中按照树形结构的运算逻辑从底层开始往上运算，第3行控制整个迭代的次数；第4行确定当前层数的元素对大小；第5行确定该层的螺旋因子；第8行第13行是主要的运算逻辑，称为“蝶形运算”，在这里完成被分割为左右两部分的元素组的DFT运算。以上就是FFT的大体叙述，FFT在计算大整数乘法中是一个比较重要的部分，它为大整数乘法提供了优化方案，是目前实现大整数乘法的一种优秀算法。接下来我们要讨论在C+下，利用FFT实现大整数乘法要考虑的细节问题。位元反转（Bit reversal）我们回顾一下之前在非递归形式的FFT中提到的位元反转问题，1中并没有给出位元反转的伪代码，但有提到我们可以实现时间复杂度为O(log n)的位元反转算法。更进一步，如果我们采用表映射的方式，我们还能以近乎O(1)的时间复杂度实现位元反转，只是我们可能需要比较多的内存空间，具体需要多少空间取决于计算机硬件和编译器。进行位元反转的简单算法可以用伪代码描述如下：这种通过移位来完成位元反转的时间复杂度为O(k), k为整数的位数。假如我们采用这种方法进行位元反转的话，那么在上面的FFT的非递归算法中BIT-REVERSE-COPY的时间复杂度就是O(nk)。一般来说，在大整数乘法中，要进行位元反转的元素个数是很多的，因此n远大于k，可以将k看成一个对效率影响不大的常数，BIT-REVERSE-COPY的时间复杂度可以认为是O(n).我们在上面做的只是按照传入的数值进行一次位元反转，假如要进行反转的数有n个，那么就要调用n次上面的方法，不过我们在实现时可以将上述算法改写成内联函数以提高效率。如果不采用上面那种逐个计算的方法的话，还有下面这种方法。其中a是长度为length的数组，算法执行完后，从a0到alengt-1都保存了与下标对应的位反转的值。这种方法计算的次数会减少很多，但是需要O(n)的空间来保存计算出的数据。我们可以通过直接读取a数组中的数据来获得位反转的值，这就是表映射的基本思想。因为要计算位元反转的元素的个数很多，假如我们能提高位元反转算法的效率的话，将会带来一定程度的性能提升。高效的位元反转算法要综合考虑计算机硬件和编译器的优化策略，在进行编程时我们可以根据具体情况，对位元反转进行相关的优化。一般来说，采用表映射的方法来实现位元反转已经能够满足性能需求了，在下面的算法实现中我们就采用了这种方案。如果要追求极致的性能的话，我们还可以进行C+代码内嵌汇编，根据CPU的指令集进行更进一步的优化。相关问题的讨论已经超出了本文的范围，如果你有兴趣，可以参考8.C+算法实现以下是C+算法实现的清单，代码在Visual Studio Community 2015下编译通过。#define _USE_MATH_DEFINES#include using namespace std;/ 复数类class Complexpublic:/ 构造器Complex(long double real = 0.0, long double imaginary = 0.0) this-real = real;this-imaginary = imaginary;/ 重载运算符inline Complex operator +(const Complex &anotherComplex) return Complex(this-real + anotherComplex.real, this-imaginary + anotherComplex.imaginary);inline Complex operator -(const Complex &anotherComplex) return Complex(this-real - anotherComplex.real, this-imaginary - anotherComplex.imaginary);inline Complex operator *(const Complex &anotherComplex) return Complex(this-real*anotherComplex.real) - (this-imaginary*anotherComplex.imaginary),(this-real*anotherComplex.imaginary) + (this-imaginary*anotherComplex.real);/ 为了操作方便和提高效率，就不把它们声明为私有了long double real, imaginary;/ 这些变量放在外面是为了方便，编译器会自动帮我们分配堆空间const int MAX_LENGTH = ;char str1MAX_LENGTH / 2, str2MAX_LENGTH / 2;Complex aPointValuesMAX_LENGTH, bPointValuesMAX_LENGTH;long long sumMAX_LENGTH;unsigned long long *bitReverseLookupTable;/ 以2的幂为基准进行对齐size_t alignByPowerOf2(const size_t lengthOfStr1, const size_t lengthOfStr2) size_t roof = 1;while (roof (lengthOfStr1 1) | (roof (lengthOfStr2 1) roof = 1;for (size_t i = 0; i lengthOfStr1; i+) aPointValuesi = Complex(str1i - 0, 0);for (size_t i = 0; i lengthOfStr2; i+) bPointValuesi = Complex(str2i - 0, 0);return roof;/ 点值表达式乘法运算void pointValueMultiply(size_t length) for (size_t i = 0; i 1;for (int k = 2; k 1 1;bitReverseLookupTablek + 1 = bitReverseLookupTablek bitReverseLookupTable1;/ 位元反转交换void bitReverseShuffle(Complex values, size_t length) for (size_t i = 0; i length; i+) / 避免元素自身交换和再次交换if (i bitReverseLookupTablei) swap(valuesi, valuesbitReverseLookupTablei);/ 迭代形式的FFT，基本上按照算法导论的来做void FFT(Complex sampleValues, size_t length, bool inverse) bitReverseShuffle(sampleValues, length);for (int i = 2; i = length; i = 1) / DFT和反向DFT中e的幂不同，根据欧拉公式需要进行正负转换long double tempValue = (inverse) ? (1) : (-1) * 2 * M_PI / i;/ 根据欧拉公式做出的变换Complex wn(cos(tempValue), sin(tempValue);for (size_t j = 0; j length; j += i) Complex w(1, 0);/ Butterfly operationfor (size_t k = j; k 1); k+) Complex u = sampleValuesk;Complex t = w*sampleValuesk + (i 1);sampleValuesk = u + t;sampleValuesk + (i 1) = u - t;w = w*wn;/ 反向DFTif (inverse) for (int i = 0; i length; i+) sampleValuesi.real /= length;void adjustResult(size_t size) / 取整for (size_t i = 0; i 0; i-) sumi - 1 += sumi / 10;sumi %= 10;void outputResult(size_t resultSize) / 关闭流同步以提升输出效率ios:sync_with_stdio(false);for (size_t i = 0; i resultSize; i+) cout sumi;cout str1 str2).fail();/ 清理点值表达式数组void flushPointValues() memset(aPointValues, 0, MAX_LENGTH*sizeof(Complex);memset(bPointValues, 0, MAX_LENGTH*sizeof(Complex);int main() size_t length;size_t lengthOfStr1, lengthOfStr2;/ 输入数据while (inputString()lengthOfStr1 = strlen(str1), lengthOfStr2 = strlen(str2);length = alignByPowerOf2(lengthOfStr1, lengthOfStr2);/ 申请位元反转映射表的存储空间bitReverseLookupTable = new unsigned long longlength;/ 计算位元反转映射表computeBitReverseLookupTable(length);/ 利用FFT算法进行DFTFFT(aPointValues, length, false);FFT(bPointValues, length, false);/ 点值乘法运算pointValueMultiply(length);/ 利用FFT算法进行反向DFTFFT(aPointValues, length, true);/ 调整结果adjustResult(length);/ 输出结果outputResult(lengthOfStr1 + lengthOfStr2 - 1);/ 清理flushPointValues();参考资料：1 Thomas H. Cormen, Charles E. Leiserson, Ronal