FFT详解&大数乘法

分享一个写的非常好的博客 http://blog.csdn.net/u013351484/article/details/48739415 http://blog.csdn.net/u013351484/article/details/48809943 比起一些写的看都看不懂的文章好到不知到哪去了

传统的乘法的方法类似于利用列竖式的方法，时间复杂度为 O(N2) 。但是利用FFT的方法，我们可以把时间复杂度降到 O(NlogN) 。

设A和B是两个很大的数，C=A*B。把这两个数转换为类似幂级数的形式，即：（不足位数的补足前导0） A=a0+a1x+a2x2+……+an−1xn−1 B=b0+b1x+b2x2+……+bn−1xn−1 （接下来所说的都以A为例，B同理） 其中 a0、a1...an−1 分别表示A的第一位、第二位…第n位。这样我们就能用一个向量 (a0,a1,...,an−1) 来表示A，这种表示方法叫做系数表示法，当x取10时就是我们常用的表示法。

我们发现这个幂级数的x其实是可以任取的，令： A(x)=a0+a1x+a2x2+……+an−1xn−1 对于这个n-1次的函数，最多只需n个点即可唯一确定这个函数，也就是说我们取这n个点： {(0,A(0)),(1,A(1)),...,(n−1,A(n−1))} 即可唯一确定这个函数了，继而确定这个数A，换句话说A这个数唯一地对应了这n个点的集合，根据这n个点的集合也能反推出A。这种表示方法叫点值表示法。

这样我们就找到了一种新的计算方式，令A、B分别由： {(0,A(0)),(1,A(1)),...,(n−1,A(n−1))} {(0,B(0)),(1,B(1)),...,(n−1,B(n−1))} 得到，将对应点的y坐标相乘： {(0,A(0)∗B(0)),(1,A(1)∗B(1)),...,(n−1,A(n−1)∗B(n−1))} 显然这就是C的点值表达式。考虑到两个n位数相乘的结果长度会加倍，所以一般都要取2n个点进行计算。

下面分析这个算法的时间复杂度： 1、得到两个数的点值表达式，因为每个点都要 O(N) 次计算，所以这一步的复杂度是 O(N2) ； 2、点值对应相乘，只需做 O(N) 次； 3、还原成系数表达式，一般来说只能采用高斯消元法，时间复杂度还是 O(N2) 。 所以我们是用更高级的方法得到了一样复杂度的做法吗？ 其实不是的，这个算法是根本无法实现的。。。 先不说还要用高斯消元法，单单是得到n个点就无法实现：在整数集合内根本无法接受形如 xn−1 的计算还不能取模。 但是，点值表达法却是FFT的关键一步。

所有有关DFT和FFT所需的约定都会放在这里，留心一下即可。 1、A和B是两个乘数，C是乘积。A和B的位数都必须是2的幂次，如果不是2的幂次就补足前导0，这是为了保证之后的二分能够正确进行。又因为C的位数要加倍，所以A和B的位数还需加倍，设这个位数为n。

在复平面上任何一个复数z都能表示成为一个向量，即： z=r(cosθ+isinθ) 其中r是z的模长， θ 是向量与x轴的夹角，称之为幅角。 定义： eiθ=cosθ+isinθ 则有： z=r∗eiθ 由此可以推出： (cosθ+isinθ)α=(cosαθ+isinαθ) 这就是棣莫弗公式。

在复数集下满足方程 xn=1 的解一共有n个，这n个解构成1的n次单位根。并且这n个根中存在至少一个根 wn 使得 wn 的1~n次恰好就是这些n次单位根，称 wn 为本原根。 换句话说： w0n、w1n、w2n、...、wn−1n 这n个数互不相同，并且这n个数的n次方都是1。 如果我们把这n个复数在复平面上表示出来，我们发现这n个复数恰好平分360度。这个东西可以用棣莫弗公式来证明，因为复数的1次方、2次方、3次方···恰好就是幅角的1倍、2倍、3倍··· 由此我们可以得到一个通用的本原根：（说通用是因为有些情况下本原根不止一个） wn=cos2πn+isin2πn

之前提到x的选取是任意的，而DFT的神奇之处在于x取的就是n次单位根。然后我们能把A写成这样： A(x)=∑n−1j=0aj∗xj 之前已经说了要取n个点，所以令： yk=A(wkn)=∑n−1j=0ajwkjn （latex不让打两次指数，相当于本原根的k次的j次） 其中 0y0,y1,y2,...,yn−1} （点值表示法）的转换。 称y为a的离散傅里叶变换(Discrete Fourier Transform,DFT)。

遗憾的是，虽然DFT看起来更加高级，但依然改变不了它 O(N2) 的事实，事实上它连 O(N2) 都不如，因为大量的浮点数运算会带来更大的常数。 接下来的FFT就能在 O(NlogN) 的时间里完成从a到y的变换。

构造两个全新的次数界为n/2的多项式： A[0](x)=a0+a2x+a4x2+...+an−2xn/2−1 A[1](x)=a1+a3x+a5x2+...+an−1xn/2−1 显然有： A(k)=A[0](k2)+k∗A[1](k2) 这样，原问题： A(x) 在 w0n、w1n、w2n、...、wn−1n 就转化为： A[0] 和 A[1] 在 (w0n)2、(w1n)2、(w2n)2、...、(wn−1n)2 类似于归并排序，这样不断分裂下去直到结束。

那么现在已经的两个多项式如何递归下去？因为转化之后的问题和原问题并不是完全符合的（多了一个平方）。 然而这个问题只需套用一下之前的棣莫弗公式，复数的平方可以看做幅角加倍，而幅角加倍的结果就是复数 wn 变为了 wn/2 。 一般地： wdkdn=wkn 称为相消定理。

所以： (w0n)2、(w1n)2、(w2n)2、...、(wn−1n)2 就能转化为： w0n/2、w1n/2、w2n/2、...、wn−1n/2 这样一来只要替换本原根，子问题就和原问题一模一样啦！

现在我们已经得到了两个子问题的解，如何合并成原问题的解？ 举一个n=4的例子： y0=A(w04)=A[0]((w04)2)+w04A[1]((w04)2) y1=A(w14)=A[0]((w14)2)+w14A[1]((w14)2) y2=A(w24)=A[0]((w24)2)+w24A[1]((w24)2) y3=A(w34)=A[0]((w34)2)+w34A[1]((w34)2) 根据相消定理可以得到： y0=A(w04)=A[0](w02)+w04A[1](w02) y1=A(w14)=A[0](w12)+w14A[1](w12) y2=A(w24)=A[0](w22)+w24A[1](w22) y3=A(w34)=A[0](w32)+w34A[1](w32) 套用一下两个公式： 1、 wk+n/2n=wkn∗wn/2n=−wkn 2、 wk+nn=wkn∗wnn=wkn 就可以得到： y0=A(w04)=A[0](w02)+w04A[1](w02) y1=A(w14)=A[0](w12)+w14A[1](w12) y2=A(w24)=A[0](w02)−w04A[1](w02) y3=A(w34)=A[0](w12)−w14A[1](w12) 由此我们可以在 O(NlogN) 的时间里完成从a到y的转换，这个过程就是快速傅里叶变换。

在经过FFT之后我们可以在 O(N) 的时间里求出C的点值表达式。但是问题来了，怎样从点值表达式回到系数表达式？

构造一个范德蒙德矩阵 V 满足： ⎡⎣⎢⎢⎢⎢⎢⎢y0y1y2:yn−1⎤⎦⎥⎥⎥⎥⎥⎥= ⎡⎣⎢⎢⎢⎢⎢⎢⎢111:11wnw2n:wn−1n1w2nw4n:w2(n−1)n.........:...1wn−1nw2(n−1)n:w(n−1)(n−1)n⎤⎦⎥⎥⎥⎥⎥⎥⎥ ⎡⎣⎢⎢⎢⎢⎢⎢a0a1a2:an−1⎤⎦⎥⎥⎥⎥⎥⎥ 不妨记做 a=Vy ，而之前的FFT就可以看作是计算一次 V 。 从代数的角度来看，FFT的逆运算就是V的逆运算，而 V 的逆运算是什么？显然是V−1，也就是逆矩阵！因为有： aV−1=VV−1y ，即： y=aV−1 具体的证明就不证了，那个博客里有，扔结论： 对于第j行第k列，01],s2[N>>1]; double rea[N],ina[N],reb[N],inb[N],ret[N],intt[N]; int i,len1,len2,lent,lenres,len; int res[N>>1]; void FFT(double *reA,double *inA,int n,int flag) { if (n == 1) return; int k,u,i; double reWm = cos(2*pi/n) , inWm = sin(2*pi/n);//本原根 if (flag) inWm = -inWm; double reW = 1.0 , inW = 0.0; for (k = 1,u = 0;k < n; k += 2,u++)//奇数项和偶数项分开 {ret[u] = reA[k]; intt[u] = inA[k];} for (k = 2;k < n; k += 2) {reA[k/2] = reA[k]; inA[k/2] = inA[k];} for (k = u,i = 0;k < n && i < u; k++,i++) {reA[k] = ret[i]; inA[k] = intt[i];} FFT(reA,inA,n/2,flag); FFT(reA+n/2,inA+n/2,n/2,flag); fo(k,0,n/2-1)//合并 { int tag = n / 2 + k; double reT = reW * reA[tag] - inW * inA[tag]; double inT = reW * inA[tag] + inW * reA[tag]; double reU = reA[k] , inU = inA[k]; reA[k] = reU + reT; inA[k] = inU + inT; reA[tag] = reU - reT; inA[tag] = inU - inT; double reWt = reW * reWm - inW * inWm; double inWt = reW * inWm + inW * reWm; reW = reWt; inW = inWt; } } int main() { while (~scanf("%s%s",s1,s2)) { memset(res, 0 , sizeof(res)); memset(rea, 0 , sizeof(rea)); memset(ina, 0 , sizeof(ina)); memset(reb, 0 , sizeof(reb)); memset(inb, 0 , sizeof(inb)); len1 = strlen(s1); len2 = strlen(s2); lent = (len1 > len2 ? len1 : len2); len = 1; while (len < lent) len >1]; double rea[N],ina[N],reb[N],inb[N],ret[N],intt[N]; int i,len1,len2,lent,lenres,len; int res[N>>1]; void Swap(double &x,double &y) {double t = x; x = y; y = t;} int rev(int x,int len) { int ans = 0,i; fo(i,1,len) {ans=1;} return ans; } void FFT(double *reA,double *inA,int n,int flag) { int s,i,j,k; int lgn = log((double)n) / log((double)2); fo(i,0,n-1)//数组重排 { j = rev(i,lgn); if (j > i) {Swap(reA[i],reA[j]); Swap(inA[i],inA[j]);} } fo(s,1,lgn) { int m = (1 len2 ? len1 : len2); len = 1; while (len < lent) len