基于FFT的大整数乘法
从另一个角度理解FFT在多项式乘法中的应用
从另一个角度理解FFT（在多项式乘法中的应用)————————————————————————————————作者: ————————————————————————————————日期：ﻩ从另一个角度理解ＦFT （在多项式乘法中的应用）n ｅg ｉiz ｈa ｏ——我一个信号处理算法,怎么用来算多项式乘法了呢？很多同学学习FFT 时，看到那“一堆符号”“一堆算式”就决定“C ｔrl+W 并去背板”。
这篇文章或许能帮助那些同学从另一个角度理解FFT 的过程。
不要被下文中的某些矩阵吓到,有些只是为了表达原理方便，在实际中我们并不需要实现这些矩阵,而是实现它的功能（如左乘这个矩阵表示什么意思）,这时应更多关注它在算法中是用来干什么的。
20１6.0３.2８更新后就用黑体小写字母ａ表示向量了。
虽说 看起来可能要更清楚,但那毕竟还是更多地用于手写。
1. 关于多项式乘法我们使用系数表示一个n 次多项式，即使用序列表示多项式()1122100--=+⋅⋅⋅+++==∑n n ni i i x a x a x a a x a x A 。
我们通常称之为多项式的系数表示。
我们也可以选取ｎ个不同的数{}1210,,,,-⋅⋅⋅n x x x x 代入多项式得到n 个值 {}1210,,,,-⋅⋅⋅n y y y y ,用n 个有序数对()()()()}{11221100,,,,,,,,--⋅⋅⋅n n y x y x y x y x 表示这个多项式，其函数图像经过了这ｎ个有序数对所表示的点。
可以证明这样所表示的多项式是唯一的。
我们通常称之为多项式的点值表示。
从系数表示得到点值表示的过程称为求值，从点值表示得到系数表示的过程称为插值。
{}1210,,,,-⋅⋅⋅n a a a a已知多项式()x A 和()x B 的系数表示,我们可以得到这两个多项式的乘积()x C 的系数表示。
为了得到()x C 中ix 项的系数，考虑当j+k=i 时，()()()i k j k k j j x x a x a x a =。
快速傅里叶变换多项式乘法
快速傅里叶变换多项式乘法快速傅里叶变换（FFT）是一种快速计算多项式乘法的算法。
在计算机科学中，多项式乘法是一个十分广泛的问题，因为它有很多应用，例如信号处理、密码学、图像处理、数据压缩等等。
FFT算法的出现解决了这个问题的时间复杂度，使得计算大型多项式乘法成为了可行的任务。
1. FFT算法的原理和步骤 FFT算法是一种基于分治思想的算法，它把多项式乘法分解成两个较小的多项式乘法，然后以递归的方式继续分解直到小到可以直接计算的程度。
FFT算法的主要步骤如下：Step 1：将两个多项式的系数表达式分别展开，组成两个系数向量，由于向量在FFT算法中的操作更加方便，因此将分别展开出的系数给予两个向量。
Step 2：给向量补齐空缺值，由于$n$（两多项式次数之和）为$2$的幂，于是补齐至$2^n$个。
Step 3：对向量进行傅里叶变换，得到两个傅里叶变换向量。
这一步可以利用DFT（离散傅里叶变换）或IDFT （离散傅里叶逆变换）实现。
具体实现方式后续详细介绍。
Step 4：将两个傅里叶变换向量逐位相乘（注意：乘法运算是复数乘法，不是单纯的数乘），得到一个新的傅里叶变换向量。
Step 5：对新的傅里叶变换向量进行傅里叶逆变换（IDFT），得到的结果就是最后的多项式系数向量。
2. DFT的计算及优化思路DFT（离散傅里叶变换）是FFT算法中的重要计算基础，因此，如何快速准确地计算离散傅里叶变换的值是至关重要的。
对于长度为$n$的一个序列$a$，它的DFT计算公式如下：$$A(k)=\sum _{j=0}^{n-1}a(j)\cdot e^{-i2\pikj/n}$$其中$i$表示虚数单位，$e$表示自然常数，$k$表示频率，$j$表示时间，$n$表示序列长度。
根据公式，可以采用暴力枚举的方式计算出$A(k)$的值，但是时间复杂度达到$O(n^2)$，很不适合计算大量数据。
于是，FFT算法采用了一种基于蝴蝶运算的DFT优化方式，能够将时间复杂度降到$O(nlog(n))$：如果序列的长度为偶数，可以将序列分成两个序列，再分别进行DFT；如果序列长度是奇数，则可以将序列复制一份，增加一个零元素，然后将序列分成两部分，分别进行DFT的计算。
FFT实现大整数乘法
目录FFT实现高精度乘法 (2)绪论 (2)DFT（离散傅里叶变换） (4)Cooley-Tukey算法实现FFT（快速傅里叶变换） (4)位元反转（Bit reversal） (7)C++算法实现 (8)参考资料： (12)FFT实现高精度乘法绪论对于这个课题，我们先从多项式开始一步一步进行分析。
一个度数为n的多项式A(x)可以表示为：一个多项式可以有不同的表示方式，一种是我们常见的系数表示，另一种是点值表示。
例如上面的多项式就是以系数表示形式定义的。
它展开之后可以写成下面这样：这种方式的优点是很直观，也是我们最常用的表示多项式的方法，但是采用这种方式表示的多项式在进行大数值的乘法运算的时候，如果采用直接乘法方式，假设两个多项式的度数都为n, 那么乘法计算的时间复杂度为O(n2).如果采用另一种不那么直观的点值表示法，对于多项式的度数比较大的情况，采用点值表示法进行乘法运算可以获得较大的性能提升。
对于上面提到的多项式A(x), 它的点值表示法大概像下面这样子：其中x0~x n-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是要一一对应的，否则不能直接进行乘法运算。
利用快速傅里叶变换计算多项式乘法
A(
1 2 n1
), B(
1 2 n1
2 ( n ) ) 对所有点求值时间:
2n 2n A(2n1 ), B(2n1 )
需要降为: (n log2 n)
二、快速求值
a0 , a1 ,, an 适当选取2n+1个点与算法， b0 , b1 ,, bn 每个点求值时间可以降为:
离散傅里叶变换dft回顾我们希望计算即在n1个n1次单位根处假设a以系数形式给出系数向量的离散傅里叶变换采用普通乘法进行计算时间复杂度为离散傅里叶变换与快速傅里叶变换fft有必要利用单位复数根的特殊性质降低时间复杂度计算离散傅里叶变换处取值转化为
利用快速傅里叶变换 计算多项式乘法
Written by szh.
{( x0 , y0 z0 ), ( x1 , y1 z1 ), ( x2 , y2 z2 ), C ( x)： , ( xn , yn zn )}
时间复杂度为 ( n)
• 点值表达下n次多项式的乘法
点值表达为：
A( x)： {( x0 , y0 ), ( x1 , y1 ), ( x2 , y2 ),, ( xn , yn )} B( x)： {( x0 , z0 ), ( x1 , z1 ), ( x2 , z2 ),, ( xn , zn )}
可以得到C(x)中n+1对点值：
而C(x)是2n次多项式
表示C(x)
需要将A(x),B(x)点值表达扩展到2n+1对
• 点值表达下n次多项式的乘法
扩展点值表达为：
A( x){( ：x0 , y0 ), ( x1 , y1 ), ( x2 , y2 ),, ( x2n , y2n )} B( x)： {( x0 , z0 ), ( x1 , z1 ), ( x2 , z2 ),, ( x2n , z2n )}
FFT实现大整数乘法
目录FFT实现高精度乘法 (2)绪论 (2)DFT（离散傅里叶变换） (4)Cooley-Tukey算法实现FFT（快速傅里叶变换） (4)位元反转（Bit reversal） (7)C++算法实现 (8)参考资料： (12)FFT实现高精度乘法绪论对于这个课题，我们先从多项式开始一步一步进行分析。
一个度数为n的多项式A(x)可以表示为：一个多项式可以有不同的表示方式，一种是我们常见的系数表示，另一种是点值表示。
例如上面的多项式就是以系数表示形式定义的。
它展开之后可以写成下面这样：这种方式的优点是很直观，也是我们最常用的表示多项式的方法，但是采用这种方式表示的多项式在进行大数值的乘法运算的时候，如果采用直接乘法方式，假设两个多项式的度数都为n, 那么乘法计算的时间复杂度为O(n2).如果采用另一种不那么直观的点值表示法，对于多项式的度数比较大的情况，采用点值表示法进行乘法运算可以获得较大的性能提升。
对于上面提到的多项式A(x), 它的点值表示法大概像下面这样子：其中x0~x n-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是要一一对应的，否则不能直接进行乘法运算。
多项式乘法与信号fft的关系
多项式乘法与信号fft的关系多项式乘法与信号FFT的关系一、引言多项式乘法和信号FFT（快速傅里叶变换）是在数学和信号处理领域中常见的操作。
虽然它们看似不相干，但它们之间存在着紧密的关系。
本文将探讨多项式乘法与信号FFT之间的关系，并解释为什么使用FFT可以高效地计算多项式乘法。
二、多项式乘法的基本原理在代数学中，两个多项式的乘法是将它们的每一项相乘，并将结果相加得到一个新的多项式。
例如，给定两个多项式f(x) = 3x^2 + 2x + 1 和g(x) = 2x + 1，它们的乘积为h(x) = (3x^2 + 2x + 1)(2x + 1) = 6x^3 + 7x^2 + 4x + 1。
传统的多项式乘法算法需要进行大量的乘法和加法操作，其时间复杂度为O(n^2)，其中n是多项式的次数。
这在处理高次数的多项式时效率较低。
三、信号FFT的基本原理快速傅里叶变换（FFT）是一种将一个信号从时域转换到频域的算法。
它可以将信号表示成一系列复数的和，每个复数表示信号在不同频率上的分量。
FFT算法的基本思想是将信号分解成若干个频率相等的正弦和余弦函数的叠加。
通过对这些正弦和余弦函数进行离散采样，并进行快速计算，可以得到信号的频域表示。
四、多项式乘法与信号FFT的关系多项式乘法和信号FFT之间的关系在于它们都可以通过多项式的系数来表示。
具体来说，将多项式的系数看作是离散信号的采样值，那么多项式乘法可以转化为信号的卷积操作，而信号的卷积可以通过FFT高效地计算。
假设有两个多项式f(x) 和g(x)，它们的系数分别为a和b。
将a和b分别看作是信号f和g的离散采样值，通过对f和g进行FFT得到F和G的频域表示。
根据卷积定理，两个信号的卷积等于它们的频域表示的乘积。
因此，可以通过将F和G相乘得到两个多项式的乘积的频域表示。
通过对乘积的频域表示进行逆FFT，可以得到多项式乘法的结果。
五、为什么使用FFT高效计算多项式乘法使用FFT计算多项式乘法的好处在于FFT算法的时间复杂度较低。
基于超大点数FFT优化算法的研究与实现_高立宁
(1)
其中 n 0 = 0,1, , M - 1 , k 0 = 0,1, , L - 1 。 由式(1)可推导：
Z (n 0 , k ) = e
- j 2 pn0 N
´e
- j 2 pn0 N
´ ´e
- j 2 pn0 N
(2)
- j 2 pn0
由式(2)可以看出，铰链因子与 n 0 和 k 两个变量 有关， 且可以将第 k 行的铰链因子分解为 k 个 e N 相乘，因此，在存储铰链因子时可以只存储第 1 行 的铰链因子，而其它行的铰链因子可由第 1 行的铰 链因子计算得出。虽然采用计算得出铰链因子会引 入新的额外乘法运算，但该处理方式可以大大节省 内存资源，这对在资源有限的情况下实现超大点数 FFT 十分重要。为了得出新的额外乘法运算对超大
第4期
高立宁等： 基于超大点数 FFT 优化算法的研究与实现
999
理方案，其采用 4 个蝶形运算单元进行并行处理， 优化了处理效率， 其缺点是增加了 4 倍的资源消耗； 文献[6]采用 Winograd 算法将大点数 FFT 拆成小点 数进行处理，虽然在处理中引入了额外的铰链因子 相乘运算和 3 次显性转置操作，但是该算法能够有 效地降低资源的消耗，能够在现有平台上实现更大 点数的 FFT 运算。 本文针对处理器硬件资源对实现超大点数 FFT 的制约，提出一种优化的 Winograd 算法实现超大 点数 FFT 的方法。 相比传统的 Winograd 处理方法， 该实现方法通过优化铰链因子存储，改变矩阵访问 方式和限制行列划分等手段使 FFT 处理能够更好 的适应处理器的架构，使其在现有平台上可以实现 更大点数的 FFT 运算。
hmax = 4 log2 N + 5
利用快速傅里叶变换(FFT)计算多项式乘法
利用快速傅里叶变换(FFT)计算多项式乘法作者:宋振华摘要本文将讨论快速傅里叶变换(FFT),利用FFT 设计一种算法,使多项式相乘的时间复杂度降低为()n n 2log Θ,以便在计算机上高效计算多项式乘法.关键词:快速傅里叶变换、多项式乘法目录一、引言 二、算法概述 三、引理1、多项式的表示(1)系数表达形式 (2)点值表达形式 (3)插值2、利用离散傅里叶变换(DFT)与FFT 导出结果的点值表达形式(1)单位复数根(2)离散傅里叶变换(DFT)(3)通过快速傅里叶变换(FFT)计算向量y3、利用FFT 计算逆DFT,将结果的点值表达形式化为系数表达形式(1)在单位复数根处插值 (2)利用FFT 计算逆DFT四、算法具体流程五、算法的实际应用:计算大整数乘法 六、参考文献一、引言基于FFT 的离散傅里叶变换(DFT)技术,是当今信息传输、频谱分析等领域中最重要的数学工具之一.在计算机编程中,我们经常需要计算两个多项式函数的乘积.对于两个n 次多项式函数,计算其乘积最直接方法所需时间为()2n Θ.本文将讨论快速傅里叶变换(FFT),利用FFT 设计一种算法,使多项式相乘的时间复杂度降低为()n n 2log Θ,以便在计算机上高效执行.二、算法概述通过FFT 进行离散傅里叶变换,将两个多项式由系数表达形式转化为点值表达形式.将这2个点值表达形式的多项式相乘,得到结果的点值表达形式.最后利用FFT 做DFT 的逆,得到结果的系数表达形式.三、引理1、多项式的表示(1)系数表达对于次数为n 的多项式()∑==ni i i x a x A 0,其系数表达是一个由系数组成的向量()Tn a a a a ,...,,10=.考虑用系数形式表示、次数为n 的多项式)(x A 、)(x B 的乘法运算. 记∑===ni ii x c x B x A x C 20)()()(. 有∑=-=ij j i j i b a c 0.可以看出,采取逐个计算i c ),20(N i n i ∈≤≤的方式进行求解,其时间复杂度为)(2n Θ. (2)点值表达 a.定义:一个次数为n 的多项式)(x A 的点值表达就是由一个有n 个点值对组成的集合)}(,,0|),{(i i i i x A y N i n i y x =∈≤≤,使得对于n k ,...,2,1,0=,所有的k x 各不相同.b.点值表达下多项式乘法对于n 次多项式)(x A ,)(x B ,设其点值表达式分别为:)(x A :)},(,),,(),,(),,{(221100n n y x y x y x y x ,)(x B :)},(,),,(),,(),,{(,,2211,00n n y x y x y x y x ，设)()()(x B x A x C =,由于)(x C 的次数为n 2,因此必须对)(x A 和)(x B 的点值表达式进行扩展,使每个多项式都包含12+n 个点值对.给定A ,B 的扩展点值表达:)(x A :)},(,),,(),,(),,{(22221100n n y x y x y x y x ,)(x B :)},(,),,(),,(),,{(,22,2211,00n n y x y x y x y x ，.则)(x C 的点值表达形式为:)},(,),,(),,(),,{(,222,222111,000n n n y y x y y x y y x y y x ，.因此,计算)(x C 点值表达式的时间复杂度为)(n Θ.` (3)插值求值运算的逆(从一个多项式的点值表达形式确定其系数表达形式)称为插值.下面将证明,n 个点求值运算与插值运算是定义完备的互逆运算.a.多项式的点值表达可以唯一确定多项式的系数表达形式.多项式的点值表达等价于矩阵方程:⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛n n n n n nnx x x x x x x x x x x x 22222121102001111⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛n a a a a 210=⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛n y y y y 210. (3-1-3-a)记系数矩阵为V ,由Vandermonde 行列式可知,()∏≤<≤-=nj i jix x V 0.而j i n j i N j i ≠≤≤∈∀,,0,,,有j i x x ≠,因此0≠V ,即V 可逆.因此对于给定点值表达,我们能够唯一确定系数表达式,且y V a 1-=.b.对于次数为n 的多项式函数()x A ,插值算法基于如下Lagrange 公式:()()∑∏∏=≠≠--=nk kj jkk j jk x x x x y x A 0)(.(3-1-3-b)容易验证,Lagrange 公式的计算复杂度为)(2n Θ.因此,n 个点求值运算与插值运算是定义完备的互逆运算,它们将多项式的系数表达与点值表达进行相互转化.2、利用离散傅里叶变换(DFT)与FFT 导出结果的点值表达形式(1)单位复数根n 次单位复数根是满足1=n ω的复数ω.n 次单位复数根恰好有n 个:对于1,...,1,0-=n k ,这些根是nk i eπ2.由Euler 公式()()θθθsin cos i e i +=可知,这n 个单位复数根均匀地分布在以复平面原点为圆心的单位圆圆周上.ni n eπω2=称为主n 次单位根,所有其他n 次单位根都是n ω的幂次.n 个n 次单位复数根0nω,1n ω,2n ω,...,1-n n ω在乘法意义下构成一个群,该群与群>+<n n Z ,同构.由此,可以得到如下推论:a.0,0,0>≥≥∀d k n ,有knnk idkdneedndk iωωππ===22. (3-2-1-a) b.*,2N k k n ∈=∀,122-==ωωnn.(3-2-1-b)c.若*,2N k k n ∈=,},0|){(2N i n i i n ∈<≤ω=},20|){(22/N i n i i n ∈<≤ω. (3-2-1-c)对推论c 的证明:由a 可知,N k ∈∀,()2/2/2k n k n ωω=.所以对于20,n k N k <≤∈,有222222/)()(k n k n n n k n n k n n k n ωωωωωω====++.得证.d.N k N n ∈∈∀,*且k 不能被n 整除,有()∑-==100n i nk nω. (3-2-1-d)对推论d 的证明:()∑-=1n i n k nω=11)(--k n n k n ωω=11)(--kn k n n ωω=11)1(--k n k ω=0. (2)离散傅里叶变换(DFT)对于n 次多项式∑==ni i i x a x A 0)(,其系数形式为T n a a a a ),,,(10 =.(3-2-2-1) 设∑=+==ni kin i knk a A y 01)(ωω,N k n k ∈≤≤,0,(3-2-2-2)则向量T n y y y y ),,,(10 =(3-2-2-3)就是系数向量T n a a a a ),,,(10 =的离散傅里叶变换. (3)通过快速傅里叶变换(FFT)计算向量y .对于n 次多项式∑==ni i i x a x A 0)(,不妨假设N p n p ∈=+,21.对于1+n 不为2的整数次幂的情况类似,此处不再讨论.FFT 采取了分治策略,采用)(x A 中偶数下标与奇数下标的系数,分别定义两个新的2n次多项式: []21-124200)(n n x a x a x a a x A -++++= , (3-2-3-1)[]21-25311)(n n xa x a x a a x A ++++= . (3-2-3-2)注意到[])(0x A 包含A 所有偶数下标的系数,[])(1x A 包含A 中所有奇数下标的系数,于是有:[])()()(2]1[20x xA x A x A +=.(3-2-3-3)所以,求)(x A 在01+n ω,11+n ω,...,nn 1+ω处的值的问题转化为:a.求次数为2n 的多项式)(]0[x A ,)(]1[x A 在点201)(+n ω,211)(+n ω,..., 21)(n n +ω处的取值.b.根据(2-2-3-3)综合结果.根据(2-2-1-c),式(2-2-3-3)不是由1+n 个不同值组成,而是仅由21+n 个21+n 次单位复数根所组成,每个根正好出现2次.因此,我们递归地对次数为21+n 的多项式)(]0[x A ,)(]1[x A 在21+n 个21+n 次单位复数根处进行求值.这些子问题与原始问题形式相同,但规模变为一半.下面确定DFT 过程的时间复杂度.注意到除了递归调用外,每次调用需要枚举T n a a a a ),,,(10 =中所有元素,将a 划分为[]},,{200n a a a a =、[]},,,,{15311-=n a a a a a ,分别与多项式)(]0[x A ,)(]1[x A 相对应.其时间复杂度为)(n Θ.因此,对运行时间有下列递推式:)()2(2)(n nT n T Θ+=.(3-2-3-4)求解该递推式,有)log ()(2n n n T Θ=.(3-2-3-5)采取快速傅里叶变换,我们可以在)log (2n n Θ时间内,求出次数为n的多项式在1+n 次单位复数根处的值.3、利用FFT 计算逆DFT 将结果的点值表达形式化为系数表达形式(1)在单位复数根处插值根据(2-1-3-a),我们可以把DFT 写成矩阵乘积a V y n 1+=,其中n V 是一个由1+n ω适当幂次填充成的Vandermonde 矩阵.⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛n y y y y 210=⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛+++++++++2121121412112111111111n n nn n n n n n n nn n n ωωωωωωωωω⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛n a a a a 210 (3-3-1-1)易知01≠+n V .即1+n V 可逆. 下面证明11-+n V),(j i 处元素11+=-+n v ij n ij ω. (3-3-1-2)考虑111+-+n n V V ),(j i 处元素ij p :()∑∑=-+=+-++=⋅+=nk i j k n n k kjn ikn ij n n p 011111ωωω. (3-3-1-3)当i j =时,1=ij p ;当i j ≠时,根据(2-2-1-d)可知,0=ij p .满足n n n I V V =+-+111.(3-3-1-4)给定逆矩阵11-+n V ,可以求出T n a a a a ),,,(10 =.其中∑=-++=n k ki n k i y n a 0111ω.(3-3-1-5)(2)利用FFT 计算逆DFT比较(3-3-1-5)和(3-2-2-2)可知,对FFT 算法进行如下修改,即可计算出逆DFT:把a 与y 互换,用1-ω替换ω,并且将计算结果每个元素除以n .因此,我们也可以在)log (2n n Θ时间内计算出逆DFT.四、算法具体流程1、加倍多项式次数通过加入n 个系数为0的高阶项,把多项式)(x A 和)(x B 变为次数为n 2的多项式,并构造其系数表达. 2、求值通过应用)1(2+n 阶的FFT 计算出)(x A 和)(x B 长度为)1(2+n 的点值表达.这些点值表达中包含了两个多项式在)1(2+n 次单位根处的取值. 3、逐点相乘把)(x A 的值与)(x B 的值逐点相乘,可以计算出)()()(x B x A x C =的点值表达,这个表示中包含了)(x C 在每个)1(2+n 次单位根处的值. 4、插值通过对)1(2+n 个点值应用FFT,计算其逆DFT,就可以构造出多项式)(x C 的系数表达.由于1、3的时间复杂度为)(n Θ,2、4的时间复杂度为)log (2n n Θ, 因此整个算法的时间复杂度为)log (2n n Θ. 五、算法的实际应用:计算大整数乘法在密码学等领域中,经常需要进行大整数乘法运算.如果对整数q p ,逐位相乘然后相加,其时间复杂度为)log (log 1010q p ⋅Θ.当q p ,规模巨大时,这种算法将会十分低效.因此,我们采取快速傅里叶变换进行优化.设01111101010a a a a p n n n n +⋅++⋅+⋅=-- ,其中N i n i N a a i i ∈≤≤∈≤≤,0,,90(5-1)01111101010b b b b q m m m m +⋅++⋅+⋅=-- ,其中N i m i N b b i i ∈≤≤∈≤≤,0,,90(5-2)令多项式01111)(a x a x a x a x A n n n n ++++=-- ,(5-3)01111)(b x b x b x b x B m m m m ++++=-- .(5-4))()()(x B x A x C =. (5-5)注意到)10(A p =,)10(B q =. 因此)10()10()10(B A C pq ==.(5-6)将大整数相乘转化为多项式的乘法,应用快速傅里叶变换,即可得出结果.六、参考文献1、《大学数学学习指南——线性代数》(第二版),山东大学出版社,刘建亚,吴臻,秦静,史敬涛,许闻天,张天德,金辉,胡发胜,宿洁,崔玉泉,蒋晓芸编著2、《大学数学线性代数》(第二版),高等教育出版社,上海交通大学数学系线性代数课程组编著3、Introduction to Algorithms(Third Edition),ThomasH.Cormen,Charles E.Leiserson,Ronald L.Rivest,Clifford Stein4、《离散数学》,山东大学软件学院,徐秋亮,栾峻峰,卢雷,王慧,赵合计编著5、《数学分析 下册》(第二版),高等教育出版社,陈纪修,於崇华,金路编著6、《高等代数》,科学出版社,丘维声编著7、《复变函数论》(第三版),高等教育出版社,钟玉泉编著8、《快速傅里叶变换乘法的性能研究》,陕西师范大学,毛庆,李顺东。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
基于FFT 的大整数乘法
1背景
对于两个长度为n 的大整数，普通的大整数乘法的时间复杂度是()2n O ，而采用一种精心构造的分治算法，则可以将时间复杂度降低为()()585.13log 2n O n O ≈。
此处则是受到快速傅立叶变换算法的启发，提出一种新的大整数乘法，该算法基于模-p 的有限域运算，采用类似于FFT 的算法，可以将大整数乘法的时间复杂度降低为()
5.1n O ，甚至，从某种意义上说，可以达到()n n O log 。
2 基础
2.1 FFT （可以参考《算法导论》及《算法概论》相关内容）
对于两个n-1次多项式
3 基于FFT 的大整数乘法
3.1大整数的表示方法
为简便起见，这里只考虑10进制大整数，但是这样并不会失去其一般性。
对于一个10进制整数，可以将其表示为：
()01221110101010a a a a A n n n n +⨯++⨯+⨯=----
这样，就可以将一个大整数对应到一个n-1次多项式上去：
()012211a x a x a x a x A N n n n n A ++++=---- ，其中{}9,8,7,6,5,4,3,2,1,0∈i a
3.2大整数的乘法
对于两个十进制大整数A N 和B N ，设A N 各个位上的数字如下：
0121a a a a n n --
而B N 各个位上的数字如下：
0121b b b b n n --
另外记多项式
()012211a x a x a x a x A n n n n ++++=---- ，
()012211b x b x b x b x B n n n n ++++=----
于是有()10A N A =，()10B N B =。
记大整数B A C N N N ⨯=，多项式()()()x B x A x C ∙=则有：
()()()101010C B A N N N B A C =⨯=⨯=
于是已知A N 和B N ，要得到C N ，可以采用如下方法：
于是，关键的问题是找到一个有效的算法以计算()()()x B x A x C ∙=。
可惜的
是，传统的FFT 算法是基于复数上的n 次单位根n
i n n ππω2sin 2cos
+=的，因此必须要进行三角函数之类的浮点运算，从而会有一定程度的误差。
那么问题在于能否找到一个合适的方法来进行整系数多项式的乘法而没有任何的误差呢？据笔者所知，还没有一种十分合适的类似于FFT 的针对于一般的整系数多项式的乘法。
但是，大整数乘法中所要用到的多项式乘法十分特殊：即两个因子多项式中，每一个系数都不大于10. 因此下面所要找的， 就是对于这样一种特殊的整系数多项式的比较合适的乘法算法。
我们的工具是模-p 的有限域。
首先，如果我们在模-p 的有限域上能够计算出： ()()()()p x B x A x C mod '∙≡
即有()220,m o d
'-≤≤≡n i p c c i i ，又根据多项式乘法可以知道，()1818199,1,0,1,0,1,0-≤=⨯≤=∑∑∑=+-≤≤=+-≤≤=+-≤≤n b a c i k j n k j i k j n k j i k j n k j k
j i ，因此如果我们挑选一个足
够大的p 使得()181->n p ，则可知p c i <，又有：p c i <',()p c c i i mod '≡，从而可知i i c c ='.
3.3关于p 的选择
于是我们现在需要做的就是对于给定的一个数p ，设计一个类似于FFT 的算法，求出()()()()p x B x A x C mod '∙≡。
我们可以将p 的取值范围缩小一点，针对于p 为素数，且p 足够大（显然要大于2n-2）的情况来讨论。
这个任务其实就是要证明在模-p 的域上的“n 次单位根”具有与复数域上的n 次单位根类似的性质。
即证明模-p 的域上的FFT 算法是合法的。
【定理1】模-p 乘法群是p-1阶循环群。
【证明】此定理等价于同余方程()p x p mod 11≡-有p-1次本原单位根。
这个定理的证明可以参考《数论概论》中关于模-p 的本原根的内容。
在一般的近世代数教材中也会有这一内容。
■
由此，记m=p-1,于是可以找到一个次单位根ϖ。
由于p>>n,所以m=p-1>n.从而A （x ）和B （x ）都可以表示成m-1次多项式。
()0111a a a A m m +++=--ωωϖ
将此()x A 分成两部分：
()()()()()
22211122121222x xA x A x a x x a x A m k k k m k k k +=+=
∑∑-≤++-≤ 于是
()()()2221x xA x A x A +=，()()()2221x xA x A x A -=-
由上式可知，要计算()x A 和()x A -，只需要先计算()21x A 和()
22x A ，然后再做一次乘法、一次加法和一次减法即可，要计算k 个自变量的函数，需要的时间是： ()()()()12/2Θ+Θ+=k k T k T ，
从而可得()()m m O m T log =。
由上面的式子可以看出，要想进行FFT 变换，其关键在于对于n 个自变量
m ωωω,...,,2
求它们的函数的问题可以转化为求
m ωωω,...,,42
对于另外一些多项式的函数问题。
但是问题是在进行逆变换时必须要满足()p m
k k mod 01≡∑=ω，根据等比数列公
式()()p m m k k
mod 0111≡--=∑=ωωωω，这样必有1=m ω，也就是说，必须找到n 次单位
根。
根据群论中的Lagrange 定理，m|p-1，从而p-1必须精心选择，使得它包含较多的素因子2.
事实上，ω必须满足两个条件：折半引理和求和引理。
求和引理已经在上一段描述过了。
折半引理指的是i m i ωω-=+2/。
因为只有这样我们才能够实施分治策略。
网络上有一种算法，是直接找出一个符合条件的大的素数p 。
这种方法有不完美之处：当n 变得很大，超出了p 的范围时，该算法就不适用了。
当然，这种方法已经能够满足现实中的大部分需求了。
而且由于是事先计算好的模-p 的域，因而节省了大量的工作，其时间性能比较好。
我们的目的是寻找一种好的方法，使得该算法能够更容易地扩展，以适应任何长度的大整数。
关于p 为素数的情况，现在貌似也找不到一种好的方法来。
对于合数的情况，参见《基于快速傅立叶变换的大整数乘法研究》这篇论文。
虽然其定理十分优美，但是没有实用性。
因为它所选取的合数c 对n 成指数级。
即12+≥n c 。
这样在做运算的时候就失去了大整数乘法的时间优势，相当于白费力。